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equipment  for  thermal  conductivity  measurements  was 
redesigned  to  meet  the  requirements  of  this  study.  An  error 
analysis  is  used  to  identify  the  data  causing  the  greatest 
amount  of  error. 

The  influence  of  coupling  on  diurnal  subsurface 
moisture  and  heat  transfer  can  be  evaluated  by  comparing  the 
transient  non-isother mal  numerical  solutions  to  the  heat  and 
water  transport  problem  under  varying  top  boundary 
conditions.  The  results  confirm  that  coupling  does  not 
appreciably  influence  the  temperature  field  but  does  affect 
the  evaporation  and  sub surf ac e  moisture  fluxes.  Thus, 
depending  on  the  direction  of  the  temperature  and  pressure 
head  gradients,  the  evaporation  and  subsurface  water  fluxes 
may  have  been  overestimated  or  underestimated  significantly 
using  an  isothermal  model.  The  analysis  of  transport 
coefficients  indicated  that  thermally  induced  moisture 
movement  becomes  dominant  at  very  low  water  contents.  It  is 
suggested  that  for  long  term  simulations  under  extended 
drying  conditions,  the  coupled  non-isot hermal  model  should 
be  preferred  in  studying  the  long  range  seasonal  influences 
on  water  transport.  The  non-isothermal  theory  can  be  applied 
to  calculate  evaporation  from  dry  land  surfaces.  Simulation 
results  show  a  very  interesting  cyclical  evaporation 
pattern . 


The  theoretical  and  field  analysis  pointed  our 
difficulties  and  deficiencies  of  non-isothermal  models,  such 


IX 


as  the  calculation  of  flow  enhancement  factors  and  liquid 
thermal  diffusion  coefficients,  the  field  performance  of 
thermocouple  psychor meter s  and  others.  Areas  of 
hydrogeological  interest  of  such  models,  such  as  the 
subsurface  waste  heat  disposal,  thermal  influences  on  the 
migration  of  radionucleides  and  other  environmental  and 
hydrological  problems  are  pointed  out. 
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CHAPTER  1 


STATEMENT  OF  THE  PROBLEM 

In  recent  years,  there  has  been  considerable  and  varied 
research  interest  in  problems  relating  to  the  simultaneous 
transfer  of  heat  and  mass  in  flow  systems.  Although  this 
progress  is  visible  in  a  broad  range  of  scientific 
disciplines  (Whitaker,  1977;  Rossen  and  Hayakawa,  1977;  de 
Vries,  1975;  Guymon  and  Luthin,  1974;  Klute,  1973,  Harlan, 
1973;  Nielsen  et  al. ,  1972;  Keey,  1972;  Fulford,  1969  ) , 
there  is  little  evidence  that  this  work  is  making  an  impact 
on  hydrogeological  research  (Stallman,  19  64,  1  967  ).  However 
there  are  a  number  of  practical  problems  that  would  be  of 
interest  to  the  hydrogeologist  concerned  with  thermal 
influences  on  the  suosurface  flow  regime,  for  example,  the 
migration  of  radionucleides  away  from  subsurface  waste 
management  sites  in  dry  desert  regions;  hydrothermal 
convection  in  porous  media;  subsurface  waste  heat  disposal; 
evaporation  from  land  surfaces;  seasonal  fluctuations  of  the 
water  table  in  cold  environments;  natural  and  artificial 
recharge  problems;  systems  for  land  treatment  of  waste  water 
involving  a  surface  soil  saturation  followed  by  evaporative 
drying;  underground  storage  of  cryogenic  fluids  and  heat 
storage  in  aquifers;  laying  of  pipelines  in  cold  regions  and 
other  frost  problems;  and  generally  the  subsurface  hydrology 
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of  arid,  semiarid,  as  well  as  cold  regions. 

Although  problems  involving  the  simultaneous  water  and 
heat  transport  have  been  studied  extensively,  many  aspects 
are  not  yet  well  understood.  The  fact  that  different 
transfer  mechanisms  operate  simultaneously,  and  often  in 
opposite  directions,  poses  serious  difficulties  in  the 
development  of  a  suitable  mathematical  framework  and 
analytical  methodologies. 


Let  us  now  examine  some  of  the  reasons  why  the  problem 
of  the  simultaneous  movement  of  water  and  heat  is  so 
complex.  Transfer  of  liquid  occurs  mainly  as  a  result  of 
gravity,  external  pressure  forces,  capillary  forces, 
sorption  forces  and  electrical  forces  (de  Vries,  1S65; 
Nerpin  and  Globus,  1969  ).  Transfer  of  vapor  is  governed  by 
several  mechanisms,  for  example,  vapor  diffusion  under  the 
influence  of  vapor  concentration  differences  and  flow  under 
the  influence  of  external  pressure,  gravity  and  friction  at 
the  liquid-vapor  interfaces.  Heat  is  transported  by 
conduction,  radiation  and  convection  (including  sensible 
heat  transfer  by  liquid  and  vapor  movement,  as  well  as 
latent  heat  transfer  by  vapor  movement  ).  The  complicated 
internal  geometry  of  the  porous  media  adds  further 
complexities  to  the  liquid- vapor- heat  transfer.  In  addition 
to  these  difficulties,  it  should  be  noted  that  even  basic 
quantities  such  as  the  soil  water  potential  or  the  hydraulic 
and  thermal  conductivities  as  functions  of  moisture  content 
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and  temperature  are  not  easy  to  measure. 
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In  general,  the  movement  of  water  and  heat  in  porous 
media  may  be  analysed  using  three  levels  of  analysis: 
molecular,  microscopic  and  macroscopic.  For  the  molecular 
level  of  analysis,  one  develops  transport  theories  based  on 
the  movement  of  molecules.  Statistical  mechanical  concepts 
might  be  used  to  overcome  the  difficulty  of  discontinuous 
media.  At  the  microscopic  level,  one  may  utilize  the 
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their  volume  averages.  These  averages  are  taken  over  an 
elementary  volume  large  enough  to  contain  a  representative 
assortment  of  pores,  but  also  small  enough  to  maintain 
approximately  constant  values  of  macroscopic  variables 
within  the  element  (Klute,  1969  ).  Therefore,  at  this  level 
of  analysis,  overall  macroscopic  values  of  physical 


properties 

of  a 

representat ive 

vo lume 

element  are 

introduced. 

such  as 

hydraulic  and 

thermal 

conduc  ti vies , 

di ff usi vities  and  volumetric  heat  capacities,  among  others. 
A  study  of  the  connection  between  the  description  of  the 
transport  processes  on  the  microscopic  and  macroscopic 
levels  provides  further  insight  into  the  transport 
processes,  especially  with  respect  to  the  calculation  of  the 
macroscopic  transport  coefficients  from  appropriate 
microscopic  or  molecular  aspects  of  the  porous  media. 

Two  main  macroscopic  models  generally  form  the  basis 
for  predicting  water  and  temperature  distributions  in  porous 
media:  the  mechanistic  model  of  Philip  and  de  Vries  (1957  ) , 
based  on  diffusion  and  heat  flow  theory,  which  has  been 
adapted  for  porous  media;  and  the  phenomenological  approach 
based  on  irreversible  thermodynamics  (deGroot  and  Mazur, 
1962;  Prigogine,  1967;  Fitts,  1962;  Katchalsky  and  Curran, 
1965;  Hanley,  1969  ).  In  both  these  macroscopic  theories, 
the  principal  problem  is  to  express  the  microscopic 
phenomena  of  flow  in  proper  macroscopic  terms,  insuring  that 
those  mechanisms  of  primary  importance  are  included  while 
those  of  secondary  importance  are  neglected.  This  evaluation 
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of  the  relative  significance  of  the  mechanisms  tends  t 
arbitrary  depending  on  the  interests  and  objectives  o 
researcher.  Because  the  transfer  mechanisms  are  not 
completely  understood,  more  complex  models  --  whic 
presumably  more  correct  in  some  details  --  may 
necessarily  be  more  correct  overall  (Gupta  and  Churc 
1973  ) . 

The  paper  published  by  Philip  and  de  Vries  ( 
consolidating  most  of  the  previous  knowledge  on 
influence  of  temperature  gradients  on  moisture  movement 
the  groundwork  for  a  better  understanding  of  simulta 
heat  and  water  movement  in  the  soil.  The  ’’simple” 
diffusion  theory  based  on  the  modified  Pick's  equation 
porous  media,  predicted  the  vapor  flux  density  as 
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vapor  in  air  (  cm2/sec  ) 

is  a  mass-flow  factor  arising  from  differences 

in  boundary  conditions  governing  air  and  vapor 

components  (  dimensionless  ) 

is  the  density  of  water  vapor  (  gm/cm3  ) 

is  the  diffusion  coefficient  of  water  vapor  in 

the  porous  medium  (  cm 2/sec  ) 


The  transport  of  water  vapor  in  soils  under  isothermal  and 
solute-free  conditions  has  been  studied  experimentally  by 
Jackson  (1964a,  b,  c,  1965)  and  Rose  (1963a,  b)  who  found 
that  (1-1)  describes  isothermal  water  vapor  transport  in  a 
satisfactory  manner. 


However  (1-1  ),  while 
conditions,  underestimated  the 
under  temperature  gradients,  so 
observed  vapor  flux  to  the  va 
diffusion  theory  ranged  from  abou 
Vries,  1957  ).  Moreover,  the 
maximum  vapor  transfer  at  a  much 
the  observed  one. 
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flow  through  regions  of  vapor  and  liquid.  This  liquid- 
assisted  vapor  transfer  at  low  moisture  contents  was 
confirmed  experimentally  under  isothermal  conditions  by  Rose 
(1963a,  b,  1965  ).  Such  liquid-vapor  interaction,  therefore, 
is  not  restricted  to  temperature  induced  diffusion  as  Philip 
and  de  Vries  suggested. 
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this  gradient 
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than  the  average  gradient  across  the  porous  medi 
whole,  because  of  the  lower  thermal  conductivity  of  t 
According  to  Philip  and  de  Vries,  the  ratio  o 
temperature  gradient  in  the  air-filled  pores  to  the 
temperature  gradient  may  be  as  large  as  3.2  for  dry  s 
higher  (Woodside  and  Kuzraak,  1958  ).  This  represents 
increase  in  the  thermally  induced  vapor  flow  because 
non-linear  dependance  of  vapor  pressure  on  temperatur 
ratio  can  be  measured  according  to  the  method  of  de 
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(1963)  when  certain  assumptions  are  made  about  the  shape  o 
the  soil  particles  and  their  thermal  conductivities. 

Taxing  into  account  the  above  factors  in  the  vapor  flu 
equation,  Philip  and  de  Vries  developed  an  approximat 
analysis  which  generally  provides  satisfactory  agreemen 
with  the  experimental  observations.  In  order  to  describ 
thermally  induced  liquid  movement,  Philip  and  de  Vrie 
modified  Darcy's  equation  by  incorporating  the  temperatur 
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porous  media  (Groenevelt  and  Bolt,  1969;  Kay  and  Groenevelt, 
1974;  Groenevelt  and  Kay,  1974  ).  Such  an  approach  was  first 
attempted  for  porous  solids  by  Luikov  (1966)  and  Taylor  and 
Cary  (1960  ).  For  a  system  in  which  the  simultaneous  flow  of 
heat  and  water  only  are  considered,  the  appropriate  fluxes 
may  be  written  in  a  form  similar  to  those  derived  by  Philip 
and  de  Vries  (Table  1-2  ).  Jury  (1973)  presented  a 
derivation  and  a  comparison  of  such  equations  with  those  of 
Philip  and  de  Vries. 


The  non -equilibrium  thermodynamics  approach  is  more 
general  than  the  Philip  and  de  Vries's  one  because  it  does 
not  invoke  any  models  or  mechanisms  and  can  be  easily 
extended  to  include  other  influencing  factors,  such  as 
osmotic  or  concentration  gradients,  among  others.  Certain 
theorems  and  restrictions  allow  the  use  of  the  Onsager 
reciprocity  relation  which  states  that 

~  I'CjW  (1-18) 

thus  reducing  the  number  of  the  unknown  phenomenological 
coefficients  (Table  1-2)  by  one.  Cary  (1963,1964  ),  among 
others,  has  demonstrated  the  validity  of  the  Onsager 
relation  in  soils.  The  main  problem  with  this 
phenomenological  approach  is  that  it  does  not  provide 
insight  into  the  transport  coefficients,  which  makes  their 
calculation  difficult.  Other  limitations  include  the 
following:  the  domain  for  which  the  equations  apply  must  be 
very  close  to  equilibrium  for  the  Onsager  reciprocal 

be  valid.  The  condition  for  linear 


relations  to 


valid 
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phenomenological  relations  inherent  in  the  general  theory 
implies  constant  coefficients,  a  situation  which  does  not 
correspond  to  real  porous  media  flow  conditions.  Although  in 
recent  years  significant  progress  has  been  made  with  regard 
to  far -from- egui librium  thermodynamics  (Glansdorff  and 
Prigogine,  1971;  Prigogine  and  Nicolis,  1977  ),  such 
progress  has  not  yet  reached  the  area  of  flow  through  porous 
media. 


The  simi 
has  prompted 
approaches  by 
Miller,  1966 
these  scienti 
advantage  of 
observed  data 


larity  of  the  equations  in  Tables  1-1  and  1-2 
a  comparison  of  the  two  previously  discussed 
several  researchers  (Gee,  1966;  Dirksen  and 
;  Hadas,  1968;  Cassel  et  al.,  1969  ).  However, 
sts  were  not  able  to  demonstrate  any  clearcut 
one  method  over  the  other  in  approximating  the 


Analytical  solutions  of  the  non 
equations  (1-6)  and  (1-7)  or  (1-16)  and  ( 
and  1-2,  respectively,  are  not  avail 
assumed  that  the  transport  coefficients  a 
equations  are  constant.  This  assumption  i 
analytical  solutions  presented  by  Cran 
and  PecK  (1958  ) ,  Luikov  and  Michailov  ( 
al.  (1971  ).  However,  the  transport  c 
conductivities  and  d if f usivit ie s,  whic 
thermophysical  properties  of  porous 
constant  but  depend  on  moisture  conten 
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Therefore,  numerical  methods  are  the 
appropriate  means  for  solving  these  non-linear  equations. 

The  numerical  approach,  however,  is  limited  by  the  lack 
of  satisfactory  measurements  of  porous  media  properties 
under  non- is othermal  conditions,  which  are  needed  for  model 
descriptions  (Fulf ord , 1 96 9 ;  Wooding  and  Morel-Se ytoux, 1976 ; 
Whit aker, 1 977  ).  Therefore,  the  difficulty  in  comparing 
theoretical  and  model  results  with  experimental  data  is  the 
major  obstacle  to  carrying  out  numerical  calculations  on 
coupled  systems. 


Exper i me 
heat  transpo 
the  laborator 
Little  or  no 
situations.  T 
one  example 
very  pronounc 
transport  in 
ground,  seas 
temperature 
vertical  dire 


ntal  investigati 
rt  have  Deen  al 
y  (Stallman,  19 
emphasis  has  bee 
he  hot  and  dry  s 
where  temperatur 
ed  because  of 
the  vapor  phase 
onal  climatic 
gradients  whic 
ction. 


ons  of  simultaneous  water  and 
most  exclusively  carried  out  in 
67;  Jackson  et  al.  ,  19  75  ). 

n  placed  on  investigating  field 
urface  layers  of  Dare  soils  are 
e  effects  on  water  transfer  are 
the  predominance  of  moisture 
.  Also,  in  deeper  layers  of  the 
trends  may  create  persistent 
h  can  transport  water  in  a 


However,  few  investigations 
Philip  and  de  Vries  theory 
field.  Investigations  of  the 
the  upper  15  cm  of  a  bare 
that  the  thermally  induced 


based  primarily  on  the 
have  been  carried  out  in  the 
non- is  ot  he  rma  1  water  flux  of 
soil  (Pose,  1968a,  b,  c)  showed 
vapor  flux  was  of  coiparaole 
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magnitude  to  the  liquid  flux  and  was  oscillatory  in 
direction  and  magnitude.  A  field  and  modeling  study  of  heat 
transfer  in  bare  soils  (Wierenga  and  deWit,  1970)  showed 
close  agreement  between  predicted  and  observed  soil 
temperatures  for  a  wet  soil.  A  study  of  the  simultaneous 
water  and  heat  fluxes  within  the  top  10  cm  of  a  field  soil 
(Jackson  et  al.  ,  1974,  1  975)  showed  that  the  Philip  and  de 
Vries  theory  best  predicts  the  soil  water  flux  under  diurnal 
conditions  at  intermediate  water  contents. 


At  present,  a  wide  discrepancy  seems  to  exist  between 
the  'field-scale'  approach  of  the  hydrogeologists  and  the 
'local  scale'  approach  of  the  soil  physicists.  Tae 
hydrogeologists'  approach  has  been  largely  empirical  and 
crudely  averaging,  while  the  more  physically  based  "micro" 
approach  of  soil  physicists  has  so  far  not  been  successful 
in  describing  the  field  condition  and  its  associated 
heterogeneity  (S wartzendruber  and  Hillel,  1973  ).  These  two 
approaches  are  parallel  lines  of  attack  that  have  not 
converged.  If  both  approaches  moved  towards  a  common  scale 
of  system  characterization,  the  above  mentioned  discrepancy 
might  cease  to  exist. 


In  order  to  reconcile  these  approaches,  I  propose  the 
development  of  a  theory  which  will  treat  the  subsurface  flow 
domain  as  an  integrated  entity,  rather  than  as  a  series  of 
isolated  or  loosely  connected  local  systems.  Several 
scientists  (Childs,  1960;  Stallman,  1964,  1967;  Freeze, 


196Sa;  Klute, 


1969;  van  Bavel,  1969;  Vachaud  et  al.  ,  1975) 
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Hydrology,  hydrogeologists  have  tended  to  forget  that  events 
at  the  surface  and  underground  waters  are  connected  —  not 
just  as  a  matter  of  geologic  continuity  —  but  in  a  real 
sense,  in  that  there  is  continuous  movement  of  water  from 
one  domain  into  the  other.  This  connecting  link  constitutes 
the  unsaturated  zone. 


Accordingly,  the  following 
forth  for  this  thesis: 

to  develop  and  apply  a  methodo 
heat  and  water  tran 
applications; 
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to  evaluate  the  influence 

subsurface  moisture  and 
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from  dry  land  surfaces. 
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CHAPTER  2 


THEORETICAL  DEVELOPMENT 


2.1  Introduction 


m  this  chapter,  the  equations  for  water  and  heat 
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1960)  : 

Continu 

it 

y 
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3p/c)t 
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(Conser 

va 

ti 

on 

of  mass 
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(2 
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(2 
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Thermal 

e 

ne 

rg 

y  eqn. 

^>c  DT/Dt  = 

T(2p/aT)e  Vv+  J4pv 

(2 

-3) 

where 

e 

is 

t 

he  mass 

densit 

y  (gm/cm3  ) 

V 

is 

t 

he  micro 

scopic 

vec  to 

r  flow 

vel 

ocit  y 

(  cm 

/s 

ec 

) 

/ 

is 

t 

he  dynam 

ic  vis 

cosi ty 

(  poi 

se  ) 

p 

is 

t 

he  point 

fluid 

press 

ure  ( 

dyne/cm2  ) 

c 

is 

t 

he  speci 

f ic  he 

at  cap 

acit  y 

(  ca 

1/gm/C 

) 

16 
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9cd 

is  t 

he 

conduc  tive 

heat 

flux  (  ca 

is  t 

he 

viscous  dis 

si 

pa 

tion  funct 

D/Dt 

is  t 

he 

substantial 

t 

im 

e  derivati 

and  all 

oth 

er 

quantities 

ar 

e  as  def 

detailed 

deri 

vation  of  th 

e 

s 

im  ul ta  neou 

momentum 

t  ra 

ns 

fer  in  porou 

s 

me 

dia  from  b 

recently 

been 

given  by  Whit 

ak 

er 

(  1977)  .  H 

of  simplicity 

a 

nd  brevity. 

a 

less  rig 

deriving 

the 

re 

quired  equat 

io 

ns 

will  be  f 

Before 

proceeding,  i 

t 

i 

s  essenti 

paramet e 

rs  of 

i 

nterest  in  t 

hi 

s 

study  an 

made  in 

deri 

VI 

ng  appropria 

te 

e 

quatio  ns , 

for  the 

det 

er 

mination  of 

t 

hese  par 

po  tentia 

1  o 

r 

pressure  h 

ea 

d 

of  inter 

tempera t 

ure  o 

f 

the  porous  m 

ed 

iu 

m,  as  func 

time , 

are 

the  main  qua 

nt 

it 

les  sough 

develop  me  nt . 

Th 

e  pressure  p 

ot 

en 

tial  of  in 

the  frac 

t  io  n 

of 

its  energy 

pe 

r 

unit  quant 

results 

from 

a 

pressure  tha 

t 

is 

different 

pressur  e 

of 

one  standard 

at 

mo  sphere. 

taken  as 

the 

un 

it  quantity 

of 

water. 

po tenti a 

1  ha 

s 

the  unit  o 

f 

le 

ngth  (cm) 

pressur  e 

head 

4 

<  (see  also 

P- 

26).  The 

interstitial  water  is  positive  in  satu 
zero  at  the  water  table  and  negative  in 
This  negative  pressure  potential  is  also 
capillary  potential  (Rose,  1966,  Corey  et 


l/cm2/sec  ) 
ion  (  sec-2  ) 
ve 

ined  previous 
s  heat,  mass 
asic  principle 


owe ver 

for 

re 

or 

OUS 

a 

ppro 

ac 

ollowe 

id 

here 

• 

al 

to 

exam 

in 

d 

the 

assu 

mp 

wh 

ich 

ar 

e  ne 

ce 

ameter 

:s. 

P 

re 

St 

itia 

.1 

wat 

er 

ti 

ons 

of 

spa 

ce 

t 

in 

t 

he 

pr 

te 

r  sti 

.  ti 

al  w 

at 

it 

y  of 

wate 

r 

f 

r  om 

th 

e  re 

fe 

If 

a  u 

mi 

t  we 

ig 

th 

e  n 

th 

e  p 

re 

an 

d  is  o 

f  te  n 

c 

P 

ressur 

e  h 

ea 

ra 

ted 

po 

rous 

m 

un 

sa tura 

ted 

m 

k 

nown  a. 

s  ma 

tr 

a 

1.  , 

19 

67  ) 

• 

ly.  A 
an  d 
s  has 
asons 
h  in 


e  the 
tion  s 
ssar  y 
ssur  e 
and 
and 
esen  t 
er  is 
that 
re  nee 
ht  is 
ssur  e 
ailed 
d  of 
edia , 
edia. 
ic  or 
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2.2.2  Assumptions 


Th 

e 

fo 

llowing 

ar 

e 

th 

e 

ma 

jo 

r 

as 

sum 

Pt- 

ions 

an 

d 

app 

roxi 

ma 

tion 

s  under 1 

yin 

g 

the 

de 

ri 

va 

ti 

on 

of 

the  co 

uple 

d 

non 

-iso 

th 

er  ma 

1  equatio 

ns 

o 

f  w 

at 

er 

an 

d 

hea 

t 

t ransfe 

r  i 

n 

porous 

me 

dia 

to  be  pre 

sen 

te 

d  in 

t 

hi 

s 

ch 

ap 

ter ; 

1) 

the 

syst 

em  is 

CO 

mp 

rise 

d 

o 

f 

a 

non- 

def 

o  r 

mab  le 

an 

d 

c 

he 

mica 

lly  inert 

ma 

tr 

ix 

wi 

th 

li 

qu 

id 

wat 

e  r 

and 

it 

s 

V 

ap 

or  i 

n  the  voi 

ds 

am 

ong 

th 

e 

gr 

ai 

ns 

> 

2) 

liqu 

id 

fl 

ow  can 

be 

de 

scr  i 

be 

d 

by 

a 

D 

arcy 

-ty 

pe 

equa 

tion 

w 

hich  i 

mplies  -- 

am 

on 

g  ot 

he 

r 

th 

in 

gs 

— 

tha 

t 

the 

flo 

w 

i 

s 

lami 

nar; 

3) 

all 

flow 

parame t 

ers 

are 

u 

ni 

qu 

e 

f 

unct 

ion 

s  , 

so 

that 

hysteresis  is  avoided; 

4)  solute  effects  are  negligible;  no  water  composition 

variable  or  electrochemical  effects  are  considered; 

5)  the  liquid  water  and  its  vapor  are  in  a  continuous  local 

phase  equilibrium  in  the  pores; 

6)  the  macroscopic  heat  flux  equation  as  presented  by  de 

Vries  (1958)  with  negligible  heat  transfer  by  air 
convection,  radiation  and  viscous  generation  of  heat 
(2-55)  is  considered  to  be  a  valid  approximation  of 
porous  media  heat  flux; 

7)  the  soil  relative  humidity  is  independent  of  temperature 

at  constant  water  content; 

8)  air  pressure  remains  constant  and  air  effects  (entrapped 
,  solution  of  soil  air)  are  negligible. 


air 
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When  the  above  assumptions  are  applied  to  porous  media, 
they  generally  present  no  serious  limitations  (de  Vries  and 
Peck,  1958  ) -  Probaoly  the  most  serious  exception  is 
assumption  3)  above,  which  removes  hysteresis  from 
consideration.  It  is  well-known  that  relationships  such  as 
the  water  character ist ic  function  (  ij;  vs  0  )  and  to  a 
smaller  degree  the  hydraulic  conductivity  function 
(K  vs  0  or  K  vs  vj; )  exhibit  hysteretic  behavior,  especially 
for  the  near-surface  domain.  However,  as  Philip  (1955) 
pointed  out,  hysteresis  does  not  in  itself  invalidate  the 
use  of  subsurface  water  flow  equations,  as  long  as  the  water 
characteristic  relation  adopted  is  appropriate  to  the 
phenomenon  under  study.  In  order  to  sidestep  the  problem, 
scientists  considered  monotonically  drying  or  wetting  flow 
problems  and  avoided  materials  which  exhibit  pronounced 
hysteretic  behavior,  such  as  swelling  soils  and  others.  In 
other  cases  a  mean  wetting-drying  characteristic  curve  was 
shown  to  be  satisfactory  under  field  conditions  (Watson  et 
al. ,  1975  ) .  Rose  (1971)  concluded  that  variations  caused  by 
hysteresis,  though  real,  are  likely  to  be  negligible  in 
relatively  dry  soils  and  therefore  it  should  be  possible  to 
ignore  hysteresis  in  practice  without  serious  errors.  Royer 
and  Vachaud  (1975),  however,  concluded  that  hysteresis  was 
too  important  to  be  neglected  in  their  field  study. 


Although  problems  involving  hysteresis  have  been  solved 
numerically,  assuming  that  the  hysteretic  functions  were 
known  in  sufficient  detail  (Whisler  and  Klute,  1965;  Rubin, 
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1967;  Ibrahim  and  Brutsaert,  1968;  Freeze,  1969a,  1971  ), 
the  complexity  of  the  developed  theoretical  techniques  for 
handling  hysteretic  phenomena  quantitatively,  and  the 
uncertainties  associated  with  them  make  the  applicabili ty  of 


such  techniques  difficult  ( 
Miller,  1966;  Topp,  1969,  1971 
1971  ).  Considerable  research 
of  hysteresis,  so  that  the  tec 
be  simplified  and  effective 
behavior  of  field  soils  (Watso 
Also  evaluations  of  the  errors 
hysteresis  and  by  the  use 
required. 


Poulovassilis,  1962;  Topp  and 
;  Poulovassilis  and  Childs, 
is  still  required  in  the  area 
hniques  developed  recently  can 
ly  used  in  the  hysteretic 
et  al.,  1975;  Klute,  1973  ). 
introduced  both  by  neglect  of 
of  approximate  treatments  are 


In 

relatio 
vapor  a 
develop 
an  d/or 
and  hea 
pr incip 
are  mor 
in  that 
domains 
Finally 
es  timat 
coupled 


the  following  secti 
nships  between  fluxes  an 
nd  liquid,  as  well  as 
ed.  Transport  parameters 
briefly  discussed.  Subsequ 
t  flow  equations  are  deriv 
les.  As  will  be  pointed  o 
e  general  than  those  usual 
they  are  valid  for  both 
,  as  well  as  for  non- 
,  a  field- or iented  meth 
ing  the  various  transpo 
,  non-isothermal  equations 


ons,  the 
d  appropriate 
for  heat  tr 
or  coefficient 
ently  the  tra 
ed  by  invoking 
ut,  the  equati 
ly  found  in  th 
unsaturated  a 
homogeneous  p 
odology  is  d 
rt  coefficient 


constit  uti ve 
gradients  for 
ansport,  are 
s  are  defined 
nsient  water 
c  onservation 
ons  developed 
e  literature, 
nd  saturated 
orous  media, 
eveloped  for 
s  involved  in 
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2.3  Vapor  Transport 


By  making  the  following  assumptions 
1)  the  water  vapor  behaves  as  an  ideal  gas,  which  means  that 
the  vapor  density  ^ (  gm/cm3  )  is  related  to  the  vapor 


pressur 

e  pv  (cm 

H^O  )  by  t 

he  ideal 

gas 

1 

aw 

(2 

-4  )  , 

and 

=  (  M/RT  ) 

e^Py 

(2-4) 

wh 

er 

e 

M 

is  t  he 

molecula  r 

weight  of 

wat 

er  vap 

or  ( 

g 

m/mo 

le 

) 

R 

is  t  he 

uni versal 

gas  const 

ant 

(  er 

9/ 

mo  le 

/K 

) 

T 

is  t em 

perature  ( 

K  ) 

? 

'is  t  he 

liquid  wa 

ter  densit 

y  ( 

gm/c 

m3 

) 

g 

is  the 

acceler at 

ion  of  gra 

vity 

(  c 

m/ 

sec2 

) 

2) 

t 

he 

tot  al 

gas  pre 

ssure  P 

is 

cons 

ta 

nt  thr 

ough 

ou 

t  the 

porous 

medium , 

a 

mo 

di 

fied  Fi 

ckian  type 

expressio 

n  to 

de 

sc 

ribe 

vapo 

r 

flux 

( 

) 

in  por 

ous  media 

can  be  wri 

tten 

as 

f  o 

Hows 

(van 

B 

avel , 

19 

52 

9 

Rollins 

et  al.,  1 

954;  de  Vr 

ies , 

197 

5) 

• 

• 

=  -f  DaV  (Mg/RT) 

Vpv 

(2-5) 

wh 

er 

e 

V  =  p/ 

(P-Pv)  is 

the  mass 

flow 

f  ac 

to 

r  wh 

ic 

h  is 

a 

lway  s 

cl 

ose 

to 

unity;  f 

is  a 

po 

re 

geo  me 

t  r 

V 

f 

actor 

(d 

im 

en 

s ionles 

s  )  ,  and 

Da.  is 

t  h 

e  m 

olec  ul 

ar 

di 

ff 

usion 

CO 

ef 

ficient  o 

f  water  va 

por  in  air 

(  c 

m2/s 

ec 

)  • 

The  mole 

cular  diff 

usion  coef 

fici 

en  t 

Da 

of 

w 

ater 

vapor 

in 

air  is  a 

function 

of  tempera 

ture 

a  nd 

total 

g 

as  p 

re 

ssure 

an 

d 

ma 

y  be  ex 

pressed  as 

(Kruger  e 

t  al 

1 

97 

0) 

I 
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Da  =  c 

(Po/P)  (T/T  o)* 

(2-6) 

where  ; 

P0  is  the  standard 

atmosp 

aeri 

c  pr 

es 

sure  ( 

1 

atm) 

f 

P  is 

the  total  gas 

pressure  (a 

tm)  ;  T 

and 

T* 

ar 

e  abso 

lu 

te  ( 

Ke 

Ivin) 

temper 

atures; 

T0  =273.  1  5  K 

;  c=0 . 

217 

cm2/ 

se 

c  and 

n  = 

1.88 

• 

(2-6) 

is  an 

empirical  relati 

onship 

has 

ed  o 

n 

both  i 

so 

ther 

ma 

1  and 

non-is 

other  ma 1 

experimen 

tal 

data 

a 

nd 

on 

theo 

re 

tical 

consid 

era tions 

(de  Vries 

and  Kr 

uger 

,  19 

67 

)  . 

The  pore 

geometry 

factor 

f  i 

s  a 

mo 

dif ica 

ti 

on  d 

es 

igned 

to  account  for 

the  effect 

s  of  the  i 

ncre 

as 

e  of  the 

di 

ff 

usion 

path 

length 

that  resul 

ts  from 

the 

pr 

esence 

0 

f  th 

e 

solid 

matrix 

and  the 

effects  th 

at  res 

ul  t 

f  r  om 

t 

he  pre 

se 

nee 

of 

the 

liquid 

phase 

enhancing  f 

low  (see  p 

.  7) 

• 

Philip 

a 

nd  d 

e 

Vries 

(1957) 

expressed  f  as  follows: 

f 

=  oc  <£> 

fo 

r  0<  6* 

( 

2- 7a) 

f 

=  °<  (a  + 

)  fo 

r  6>0k 

( 

2- 7  b) 

where 

o(  is 

tort uosi 

ty 

(  di 

me  ns 

io 

nles  s 

)  , 

us 

ually 

approximated  by  the  const 

ant  va 

lue 

0.  66 

( 

Penman 

f 

1  940 

) ; 

<j>  is 

porosi 

ty  (  dimensionless 

)  and 

Ox  i 

s  t  h 

e 

moi stu 

re 

con 

te 

nt  at 

which 

liquid 

continuity 

first 

occ 

urs 

in 

the  por 

ous 

me 

dium. 

This  latter  quantity  is  a 

pproxi 

mated  he 

re 

by 

th 

e  m 

oi 

st  ur  e 

content  corresponding  to 

the  m 

at  ri 

c  po 

te 

n  tia  1 

of 

one 

b 

ar  of 

t  he  ch 

aract eri 

stic  curve 

(  ^  vs 

9  ). 

In 

eq 

ua tion 

(2-7 

b 

) ,  x 

is  a  f 

unction 

of  moisture 

conte 

nt. 

It  i 

s 

approx 

im 

ated 

a 

s  x  -» 

a/ax 

(Philip 

and  de  Vriesr 

1  95 

7  )  , 

where 

a 

-  t 

^  (ai r 

content)  and 

ax  =  ^"  k  » 

ten 

ds 

to 

un 

ity 

wh 

en 

1 

iquid 

continuity  no 

longer  exis 

ts  and 

t  o 

zero 

n 

ear  sa 

tu 

rati 

on 

• 
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The  pore  geometry  factor  for  a  gas  insoluble  in  water 
may  also  be  represented  as  (Currie,  1960,1961;  Rose,  1968a, 
b) 

f  =  (D/Ds)  (Ds/D0)  (2-8) 

where  D  is  the  diffusion  coefficient  in  moist  soil  of  a  gas 
insoluble  in  water  (  cm2/sec  )  ;  D0  and  Ds  are  diffusion 
coefficients  of  the  same  gas  in  air  and  dry  soil 
respectively  (  cm 2/sec  ).  The  ratio  (  D/Ds  )  accounts  for 
the  reduction  of  diffusion  by  the  presence  of  water,  and 
(  Ds /D0  )  accounts  for  the  tortuosity  or  the  presence  of 


th 

e  re 

duct io 

( 

bs  /Do 

)  acc 

so 

il  ma 

trix. 

In 

or  de 

var iabl 

es  of 

th 

er  mod 

ynamic 

us 

ed  (E 

dlef se 

as 

sumpt 

ion  t 

va 

por  a 

nd  lig 

he  vapor  flux  (2 

-5)  i  n 

terms  of 

as  T  and  9 

or  (jj 

,  tne 

expressed  in  eg 

uation 

(2-9)  is 

1943  ) .  It  is 

based 

on  the 

us  equilibrium 

exists 

be  tween 

Py  = 

ps  h 
rv 

=  Pj 

ex 

P£( 

Hg/RT) v|>  } 

(2-9) 

wh 

ere  ps 

is  t 

he  s 

at 

ura 

ted  water 

va 

por 

P 

re 

ss 

ur 

e 

(  c 

m 

H^O  ) 

an 

a  his 

the 

soi  1 

r 

ela 

tive  humidi 

ty 

(di 

me 

ns 

io 

nl 

es 

s  )  . 

0 

n  the 

ba 

sis  of 

ar  gu 

me  nt 

s 

pre 

sented  by  P 

hi 

lip 

an 

d 

de 

V 

ri 

es  ( 

19 

57  )  , 

it 

can 

be 

show 

n 

th 

at  the  rela 

ti 

ve  h 

um 

id 

it 

Y 

h 

is  m 

ai 

nly  a 

fu 

net ion 

of  m 

oist 

ure  c 

ontent,  and 

t 

hat 

th 

e 

sa 

tu 

ra 

tion 

vapor 

pr 

essur e 

Pv 

is 

a 

f 

unction  of 

te 

mper 

at 

ur 

e 

on 

iy 

.  Th 

us 

,  the 

va 

por  pressur 

e  ca 

n 

be 

expressed  a 

s 

P„  (0.  T)  =  P„S  (T)  h  (0) 


(2-10) 


Applying  the  chain  rule  of  differentiation  to  Vpv  in  (2-5) 
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and  making  use  of  (2-10  ),  the  following  relationship  is 
obtained: 

7Py=  p'  (Sh/sejve  +  h  (apj/dT)  (VT),  (2-ii) 

where  (VT)a  represents  the  average  temperature  gradient 

across  air-filled  pores  in  the  porous  medium.  This  gradient, 
as  is  mentioned  in  Chapter  1  (p.  7),  may  be  several  times 

larger  than  the  overall  macroscopic  gradient  Vt,  so  that 

(VT)a=^VT  (2-12) 

where  £>  1.  As  indicated  by  Philip  and  de  Vries  (1957)  and 
Rose  (  1968a,  b)  ,  typical  values  of  £  range  from  1.3  to  3.0 
for  most  soil  materials. 

In  order  to  replace  the  5h/^0  gradient  in  (2-11)  by  a 
pressure  head  gradient,  (2-9)  is  differentiated  with  respect 
to  0  to  give 

3h/£0  =  (Mg/RT)  h  (2-13) 

so  that  (2-11)  becomes 

Vpv  =  p*  (Mg/RT)h  V6  +  h(dpJ/dT)  ((VT)a/VT)  VT  (2-14) 

or 

Vpv  =  pvs  (Mg/RT)  liVf  +  e;  h  (dp*/dT)  Vt  ( 2- 1 5) 

Substituting  (2-15)  in  (2-5)  yields  the  vapor  flux  in  terms 
of  the  unknowns  ^  and  T 

qv/^  =  -fD<xV(Mg/RT)  2PyS  h  V^-  ^f  Da  V  (M  g/RT )  h(dpJ/dT)  VT  (2-16) 
By  defining  the  following  generalized  parameters: 

Kvyv  =  fDaV  (Mg/RT)  2Py  h  (2-17) 

where  Ku,v  (cm/sec)  is  the  vapor  conductivity;  and 

DTv  =  ^  f  Da  V  ( M  g  /R  T )  h  dp*  /dT  (2-18) 

where  DTv  (  cm2/sec/C  )  is  the  thermal  vapor  diffusivity. 
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the  vapor  flux,  can  be  more  simply  written  as 

qv/^  =  DTvVT  (2-19) 

For  the  saturated  zone,  the  vapor  flux  becomes  zero. 


2.  4  Liquid  Transport 


The  starting  point  for  the  theoretical  formulation  of 
liquid  transport  is  the  momentum  principle  --  expressed  by 
the  Na vier-S tokes  equation  (Schlichti ng ,  1968;  Le  Mehaute, 
1976)  --  which  describes  the  isothermal  flow  of  a  Newtonian 
incompressible  fluid  with  constant  viscosity 


^  DV/Dt  =  ^  g  -7p  +}n  V2V 

gravity^  pressure^  jfriction 
Inert i al  Forces 1  Applied  Forces 


(2-20) 


Because 
extreme 
zero  a 
must  be 
instead 
velocit 
conside 

1)  for 
which  i 

2)  the 


iy 

t 

t 

y 

re 

s 

vi 


t  he 

bounda 

ry  conditions  in 

porous 

media  are 

com 

plex  (the 

microscopic  flow 

velocity 

V  equals 

t  he 

surface 

of  every  grain  ) , 

the  abo 

ve  equation 

r  a 

ns 

form 

ed 

to 

a 

m 

ore  u 

sef  u 

1 

f 

or  m. 

Co 

nseq 

ue 

ntly , 

of 

deal 

in 

9 

wit 

h 

micros 

copi 

ca 

11 

y  vary 
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mean  velocity  of  the  fluid  flowing  through  a  porous  medium, 

which  is  an  empirical  approximation; 

the  Navier- St okes  equation  (2-20)  reduces  to 

0  =  pg-^p+ljVk  )V  =  ^(g-Vp/£ )  +  (y/k  )V  (2-21) 

where  V  is  the  macroscopic  (mean)  vector  flow  velocity  and  p 
is  the  mean  fluid  pressure.  Both  of  these  quantities  are  the 
outcome  of  smoothing  the  actual  microscopic  distribution  of 
velocity  and  pressure  over  a  volume  larger  than  that  of  the 
individual  pores.  Solving  (2-21)  for  V  one  obtains: 

v  =  -Ue/>)  C  -Vp/^+g  ]  =  -  (*?//)  \7$'  = 

-(k^g//)V<f>=  -KV4>  =  -KV(^+z)  (2-22) 

where  <£>'  =  g^=  g[  (p/^g)+z]  =  g  (vj>+z)  (2-23) 

$  is  the  total  (hydraulic)  potential  corresponding  to  tne 
hydraulic  head  (  cm  ) ; 

k  is  the  intrinsic  permeability  (  cm2  )  of  the  porous 
medium;  and 

z  is  the  elevation  (cm  ),  here  taken  as  positive  in  the 
upward  direction. 

(2-22)  is  the  Darcy  equation  (Hub bert , 1 956)  for  the  steady 
isothermal  solute-free  flow  of  water  in  isotropic  saturated 
porous  media  with  constant  hydraulic  conductivity  K.  By 
expressing  that  equation  in  the  following  form,  the 
macroscopic  flow  velocity  for  an  anisotropic  medium  is 
defined  as 


V  = 


=  q^/ja  =  -kv<£ 

where  K  is  the  hydraulic  conductivity  tensor. 


(  2-24) 


In  summary,  provided  that  the  inertial  terms  are 
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negligible,  as  is  generally  the  case  with  porous  media,  the 
macroscopic  mean  flow  velocity  is  directly  proportional  to 
the  macroscopic  mean  potential  gradient.  This  represents  the 
hydrodynamic  element  in  Darcy*s  equation.  Implicit  in  this 
equation  is  the  statistical  requirement  that  the  medium  must 
be  "su ff ic ie nt ly  homogeneous"  on  the  scale  of  the  averaging 
volume  --  the  so-called  "Darcy  scale"  --  for  local  mean 
quantities,  such  as  V  and  p  or  ,  to  exist  (Philip,  1970  )  . 


Darcy's  equation  can  be  applied  to  unsaturated  porous 
media  (Richards,  1931;  Childs  and  Collis-George ,  1950)  when 

the  hydraulic  conductivity  K  is  allowed  to  vary  as  a 
function  of  the  matric  potential  or  the  volumetric 

moisture  content  6.  This  change  is  due  to  the  fact  that  the 
unsaturated  hydraulic  conductivity  decreases  very  rapidly  as 
the  porous  media  water  content  decreases  from  its  saturation 
value  and  as  the  pressure  head  becomes  negative.  Thus,  for 
unsaturated  media,  (2-22)  is  expressed  as 

q  L/c  =  -  K  (vp)  =  -K(0)V$  (2-25) 

Assuming  the  coordinate  z  to  be  positive  in  a  downward 
direction,  the  hydraulic  head  becomes 

=  vjj  -  z  (2-26) 

so  that  (2-25)  is  written  as 

qe/£  =  -K(<J>)  V['J'(0.T)-z]=  -KT'f'  +  Ki  (2-27) 

where  both  K  and  vb  are  functions  of  moisture  content  and 
temperature,  and  i  is  a  unit  vector  in  the  positive  z 


direction. 
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In  order  to  generalize  Darcy's  equation  (2-22  and  2-27) 
to  take  into  account  non-isothe rmal  conditions,  (2-27)  is 
expanded  using  the  chain  rule  of  differentiation 

=  -K  70  -  K  9^/^T  VT+Ki  (2-28) 

or 

=  -D0£  V@  -Dt^  VT  +  Ki  (2-29) 

where  the  following  generalized  diff usivities  are 
introduced : 

Dgg  =  K  (5^/<^0)T —  liquid  moisture  dif fusivity (cm2/sec) 

(2-30) 

D-p|  =  K  (c^/^T)^ —  liquid  thermal  di ff  usivi t y  (cm2/sec/C) 

(2-31) 

In  employing  the  diffusivity  concepts  mentioned  above,  it 
should  be  stressed  that  although  vapor  transport  through 
porous  media  is  mainly  carried  out  by  true  diffusion,  tne 
process  of  liquid  movement  is  not  achieved  by  diffusion  but 
by  viscous  flow.  The  term  "diffusivity"  arose  in  order  to 
set  the  transport  equation  in  a  form  analogous  to  tne 
equations  of  diffusion  and  heat  conduction.  For  equations  of 
this  form,  numerous  mathematical  solutions  are  available 
(Crank, 1956;  Carslaw  and  Jaeger,  1959  ).  However,  the  main 
problem  in  applying  the  liquid  moisture  diffusivity  concept 
to  real  situations  is  that  such  a  concept  is  restricted  to 
homogeneous  unsaturated  porous  media  (Kirkham  and  Powers, 
1972  ).  Therefore,  in  this  work,  (2-28)  and  (2-29)  are 
written  as  functions  of  and  T 

q^  =  -Kv^yvjj-D-ftVT  +Ki  (2-32a) 
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where  =  K  (2-3 2b) 

Thus,  (2-32a)  is  more  general  than  (2-29)  because  it  is 
applicable  to  both  saturated  and  unsaturated  porous  media 
and  valid  for  heterogeneous  conditions. 
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3v|yaTlQ  =  (2/r)  d<r/dT=  (v^/cr)  d<r/dT=  (dln<r/dT)v|;  =  (2-34) 

where  ^  =  dln<r/dT  is  the  temperature  coefficient  of  surface 
tension,  which  is  reasonably  constant  in  the  range  of  1 0  to 
30  C.  Thus,  the  following  simple  surface  tension  model  was 
used  by  Philip  and  de  Vries  (1957)  and  numerous  subsequent 
investigators  (Gee, 1966;  Rose,  1968a,  b,  c;  Cassel  et  al. , 
1969;  Fritton  et  al.  ,  1970)  to  estimate  D-p^  : 

dT/  =^k  (2-35) 


However,  there  are  problems  with  this  simple  formulation. 
Measured  values  of  were  consistently  higher  than  the 
calculated  estimates  from  the  above  equation  (Wilkinson  and 
Klute,1962;  Cary,  1966;  Jury  and  Miller,  1974  ).  The 
explanation  for  the  differences  is  that  several  other- 
mechanisms,  besides  the  previously  mentioned  surface  tension 
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model,  also  contribute  to  liquid  flow  in  the  presence  of 
temperature  gradients  (Cary,  1966  ). 

In  order  to  circumvent  this  difficulty,  I  propose  an 
alternative  formulation  which  gives  a  better  approximation 
to  measured  data.  Invoking  assumptions  5)  and  7)  listed 

previously  in  this  chapter  (  see  also  2-9  and  2-10  )  ,  the 

following  relation  can  be  deduced  by  differentiation: 

9^/^T  =  (R/Mg)  lnh  (2-36) 

so  that 

Drl  =  K  (R/Mg)  1  lnh  |  (2-37) 

Both  the  above  described  model  (2-37)  and  the  surface 
tension  model  for  estimating  D^,  however,  suffer  from  the 
limitation  that  they  are  not  applicable  to  saturated 

conditions.  Taylor  and  Cary  (1960  ),  reported  several 
instances  of  thermally  induced  flow  in  saturated  porous 
media. 

2. 5  Total  Water  Transport 

In  order  to  obtain  the  combined  liquid  and  vapor  flux 
or  total  water  flux  qw/(>  ,  both  liquid  (2-32)  and  vapor 

(2-19)  flux  components  are  added  algebraically  as  follows: 
q*/^  =  qt/?+g,/^  =  -  V\|i  -  dtVt  +ki  (2-38) 

where 

K^=Ki^v+K^^  —  total  hydraulic  conductivity  (cm/sec)  (2-38) 

D-|-=D-pv  +D-p{  --  thermal  water  dif fusivity  (cm2/sec/C)  (2-40) 

The  total  n on-i so thermal  water  flux  (2-38)  is  governed  by 


Darcy's  equation  plus  terms  included  to  account  for 
flux  and  thermally  driven  liquid  and  vapor- 
constitutive  relation  was  termed  "the  general  Darcy  law 
Guymon  and  Luthin  (1974  ). 


2 -.6  Transient  Flow  Equations 
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(Cooper,  1966;  De  Wiest,  1966;  Eagleson , 1 970)  is 
in  (2-4  1  )  as : 

a(^>sw)/3t  +  V(^v)  =  o 

where  Sw  is  the  volumetric  degree  of  saturation 
By  combining  (2-41)  with  the  previously  ier 
equations,  the  transient  flow  equations  result.  S 
Darcy's  equation  (2-22)  into  (2-41)  yields 

3(^Sw)/^t  =V[  V(++z)  ] 

Expanding  the  lef t- hand-side  of  (2-42)  and  subs 
it  the  mathematical  definitions  of  water  and 
compressibilities  (Jacob,  1950;  De  Wiest,  196 
1966;  Cooley,  1971  ),  one  obtains  (2-43  ) 
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where  Ss  is  the  specific  storage  coefficient  (  cm-1  )  and 

F  =  30/^  +  (0/^)  %  (2-44) 
is  a  generalized  storage  coefficient  (  cm-1  ) .  It  is  assumed 
that  the  spatial  variations  in  ^  are  negligible  due  to  the 
slight  compressibility  of  water. 


(2-43)  is  the  transient  isothermal  flow  equation  which 
applies  to  both  saturated  and  unsaturated  porous  media.  For 
saturated  media,  (2-43)  reduces  to  the  elastic  storage 
eq  uation 


Ss  bty/Zt  =  VK[  V  (V|*Z)  ] 

For  unsaturated  media,  (2-43)  reduces  to  Richards' 

cw  dvp/dt  =  Vk[  v  (f+z)  1 

where  Cw  = 

is  the  specific  water  capacity  (  cm-1  ). 


(2-45) 
eq  ua  tion 
(2-46) 
( 2-47) 


By  applying 

the  continuity 

principle  to 

the 

non- iso t her ma 1 

vapor 

(2-19  ) ,  liquid 

(2-29)  and  total 

water 

flux  (2-38  )  , 

the 

following  differential  equations 

are 

obtained: 

^9v/at  =  -  V  (qv  /(> )  +E  =  V(K^7vl/+DTvVT)  +  E  (2-43) 

=  -  V(gt/^)-E  =  V(KvV(V^+D-^VT-Ki)  -E  (2-49) 

d0/at  =  -  V(g w/^)  =  V(Kv^^+DTVT-Ki)  (2-50) 

The  possibility  of  water  transfer  from  the  liquid  to  the 
vapor  phase  and  vice  versa  is  accounted  for  by  the 
source-sink  (evaporation-condensation)  terra  S  in  (2-48)  and 
(2-49)  (de  Vries,  1  958  ). 


Assuming  that  liquid  and  vapor  phases  are  always  in 
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equilibrium  --  which  is  shown  to  be  a  valid  assumption 
(Novak  and  Coulman,  1974)  --  the  following  relation  holds 

(de  Vries,  19  58)  for  any  small  region  of  porous  medium: 

0*  =  (f-0< )  ?v/^=  (<M«)  ( £/?> h  = 

exp(  (Mg/BT)^]  (2-51) 

Differentiating  the  above  relationship  with  respect  to  time 
by  means  of  the  chain  rule  and  substituting  the  result 
together  with  the  relation  Q  =  9j  +  9V  in  (2-50  )  ,  one 
obtains  the  general  non-iso thermal  transient  flow  equation 
£0/^t  =  2fy/3t  +  = 

C  exp  {Mg/ET}  ft  )  (M g/RT)  /ae  ]  2>9t/3 1 

+[  (  )  exp  (MgvjyRT}  (<j>- )  (1/$  dpJ/dT-yig^/RT?)  ]dT/dt 

=  V  (Ky^V'V  ~Ki)  (2-52) 

or 

d8/H  =  ( 1  +  A )  ^0/St+(B-G)  ^T/bt  =  V(Kvy^+DTVT-Ki)  (2-53) 
where 


A  =  (^/£)  exp[  Mg^/RT  ]  (j>~  9f  )  (Mg/RT)  d<p/dSf 

B  =  (  (£/<  )  exp[  Mg\^/RT  ]  if-fy  >  (i/?‘ ) 

G  =  l?v5/^)  exp[Mgv^/RT]  )  (Mg/HT2)  vjl 

However,  many  terms  in  the  above  equation  are  negligibly 
small  compared  to  dominant  terms  in  that  equation.  Thus  term 
A  in  the  equation  above  is  of  the  order  of  10~12  ) 

(  cm-1  )  ,  term  B  is  of  the  order  of  10-7  (^-$£) 

(  C-1  )  ,  and  term  G  is  of  the  order  of  10~14  @0  (  cm-1 

C~ 1  ).  Therefore,  for  moisture  contents  above  extreme 
dryness,  terms  A,  B,  and  G  can  be  neglected;  so  that  the 
above  equation  becomes  identical  to  (2-50  ). 


•  -  ,  $1  I  '  I 


In  order  to  express  (2-48  ) , 
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(2-49) 

and  (2-50)  in 

terms 

,  the 

concept  of 

the 

(2-44) 

is  invoked. 

Fo  r 

instance  (2-50)  becomes 

F  o)v|y^t  =  V(K^Vvjj+DTVT-Ki)  (2-54) 

The  above  relationship  constitutes  the  general, 
non-iso t her ma 1  transient  flow  equation  for  both  saturated 
and  unsaturated  porous  media. 

2.7  Heat  Transport 
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are  the  specific  heats  of  water  vapor  (at 
constant  pressure)  and  liquid  water,  respectively. 

The  separation  of  q^  into  these  component  parts  is 
discussed  by  de  Vries  (1958  ),  who  recognized  the 
approximate  character  of  (2-55).  The  approximate  nature  of 
that  relationship  stems  from  the  fact  that  while  all 
processes  of  heat  transfer  were  assumed  to  take  place 
uniformly  throughout  the  porous  medium  (notice  the  additive 
structure  of  2-55  ),  they  in  fact  do  not.  Furthermore,  it  is 
assumed  that  heat  transfer  by  air  convection  and  radiation 
is  negligible,  which  is  usually  the  case  (de  Vries,  1958  ). 

Identifying  the  heat  of  vaporization  per  unit  mass  at 
temperature  T  by 

L*  (T)  =  L*  (To  )  +cv  (T-T„  )-c*  (T-T0  )  (2-56) 

in  (2-55) ,  the  following  relation  is  obtained: 

=  -3*Vt  +L»  (T)  qv  gw  (T-T„  )  (2-57) 

or,  by  replacing  the  water  vapor  flux  qv  above  by  (2-19), 
one  obtains 

qw  =  -;WT  -^L'(T)  +C^qwAT  (2-58) 

where 

}  =  .}*  +^L'  (T)  DTv  (2-59) 

is  the  thermal  conductivity  of  a  porous  medium  including 
the  effect  of  heat  transfer  by  vapor  movement  under 
the  influence  of  temperature  gradients. 

Experimentally,  it  is  this  quantity  which  is 
measured  and  not  because  when  a  temperature 
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M  =  (RT/Jh)  ah/aT 


(2 -63f ) 


is  the  heat  content  of  the  dry  porous  medium;  the 

latent  heat  of  the  vapor;  Hv  and  the  heat  contents  of 

vapor  and  liquid  water,  respectively;  and  Hw  is  the  heat  of 
wetting.  ( 2  —6  3f ) —  representing  the  differential  heat  of 
wetting  W  --  is  obtained  according  to  thermodynamic  theory 
(Edlefsen  and  Anderson,  1943)  ;  J  in  that  equation  is  the 
mechanical  equivalent  of  heat.  By  grouping  the  heat  capacity 
terms  from  (2-63a,  c  and  d  ) ,  the  following  total  heat 
capacity  of  a  moist  porous  medium  may  be  defined 

C  =  c^  +cv  ^  9V  +c^  ^  9i  (2-64) 

Combining  (2-55)  and  (2-63)  in  the  heat  conservation 
equation  (2-62)  and  making  use  of  definitions  (2-56)  and 
(2-64  ) ,  one  obtains 

C  dT/3t  +  L  '  (T)  (  <>  Vqy)  -£W  ^t/^t  = 

VQ*VT)  ~cvqvVT  -ciqAVT  (2-65) 

where  the  continuity  equation,  as  indicated  by  the  first  two 
terms  of  (2-50)  is  used  as  well.  By  combining  the  last  two 
terms  of  the  right-hand  side  of  (2-65),  the  following 
general  transient  heat  flow  equation  is  obtained 


C  dT/^t  +L'  (T)  (  ^  ^0y/^t+  Vqv  )  -^W  'bQt/dt  = 

V  Q*V T)  ~cA  qwVT  (2-66) 

The  term  <jW  representing  the  heat  content  change  due 

to  liberation  of  the  heat  of  wetting  is  a  small  quantity, 
especially  for  porous  media  with  a  small  clay  fraction  and, 
therefore,  can  be  neglected.  The  ^)0v/^t  term  may  also  be 
neglected  due  to  its  negligibly  small  magnitude.  Finally 


■ 
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replacing  gv  in  the  above  equation  by  (2-19  ),  one  obtains 
the  reduced  transient  heat  flow  equation 
C  ai/^t  =  VQV T)  +(3L'  (T)  V(K^V^)  -c£qwVT  (2-67) 

At  saturation,  the  second  term  of  the  right-hand  side  of  the 
above  equation  becomes  zero  and  (2-67)  reduces  to 

C  3T/£t  =  V$VT)  -  C^qwVT  (2-68) 

Table  2-1  summarizes  the  final  forms  of  the  developed 
equations  and  transport  coefficients. 

The  derived  flow  equations  (2-54)  and  (2-67)  form  a 
system  of  coupled  partial  differential  equations  which,  in 
principle,  may  be  solved  in  order  to  predict  the  evolution 
of  a  given  flow  system  subject  to  appropriate  boundary  and 
initial  conditions.  However,  the  flow  equations  are 
non-linear;  that  is,  the  various  transport  parameters 
appearing  in  (2-54)  and  (2-67)  are  functions  of  the 
dependent  variables  vj;  (or  0  )  and  T;  this  fact  makes  the 
solution  of  the  above  equations  very  difficult.  Chapter  4 
contains  a  more  complete  elaboration  of  this  problem  and  a 
discussion  of  methods  that  may  be  utilized  to  solve  the  flow 
eq  uations. 

2^8  Field-or ien t ed  Methodology  for  Transport 

Coef f ic ie  nts 


As  mentioned  in  Chapter  1,  the  coupled  non-isoth er mal 
flow  equations  have  not  been  adequately  tested  by 
application  to  field  problems.  Accordingly,  field  data 
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necessary  to  develop  relevant  models  --  are  extremely 
scarce.  One  of  the  major  difficulties  in  applications  of 
such  models  to  field  situations  is  the  calculation  of 
transport  coefficients  because  they  constitute 

dif f icu lt-t o- mea sure  physical  properties  cf  the  porous  and 
fluid  media.  In  this  section  I  will  develop  simplified 
procedures  for  estimating  the  transport  parameters  for  a 
unified  unsaturatea-sa turated  profile  without  sacrificing 
accuracy  in  relation  to  the  analytical  relations  of  Table 
2-1.  Next,  the  heat  and  water  transport  equations  will  be 
considered  in  one-dimensional  form.  Formulations  based  on 
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literature  data  or  simply  as  indicated  in  Chapter  4  (p. 
141),  the  evaluation  of  D*  may  not  be  difficult  in  the 
field,  provided  that  the  soil  psychrome te rs  used  for 
measuring  soil  relative  humidity  (h)  function  properly. 
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diffusivity  D-j^  may  be  calculated  from  the  non-iso th ermal 
liquid  flux  equation  (2-32)  as 
Du  =  (-K^J/^Z  -qt/^  +Ki)/(ST/^z)  (2-73) 

where  the  saturated  K  may  be  measured  by  piezometric  methods 
(Chapter  3  )  . 


This  proposed  methodology  may  be  used  to  establish  and 
solve  non-isothermal  equations  for  field  systems.  Although  a 
discussion  on  the  accuracy  and  limitations  of  this  approach 
is  presented  later  (Chapter  3  ) ,  it  is  important  to  remember 
that  this  methodology  depends  mainly  on  the  accuracy  of  tne 
estimation  of  the  water  fluxes  and,  particularly,  of  the 
heat  fluxes. 


CHAPTER  3 


FIELD  AND  LABORATORY  STUDIES 


3.  1  Intro duct  ion 

A  general  lack:  of  adequate  data 
heat  and  water  transfer  in  porous  med 
field-tes ti ng  of  the  non-is ot herm al  t 
coupled  systems.  For  this  reason, 
thermally  induced  flow,  especially  un 
has  not  yet  been  definitely  assessed, 
the  set  of  field  experiments  for  th 
there  were  no  comprehensive  f 
non-isother ma 1  flow  reported  in  the 
aspects  of  the  problem  had  been  consi 
worked  with  non-i sother mal  water 
Wierenga  and  de  Wit  (1970)  worked  wit 
In  recent  years,  there  has  been  evid 
treatment  of  the  coupled  aspects  of  t 
al.  (1974,  1975)  investigated  heat 
field  conditions.  However,  their  work 
the  top  10  cm  of  the  soil,  and  seve 
in  the  non-isother mal  flow  eguatio 
measured  in  their  experiments. 
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Thus,  an  important  object 
design  a  field  experiment  to  meas 
physical  quantities  that  would 
methodology  for  non-isother ma 
conditions  that  I  proposed  in  Ch 
program  of  field  studies  will  be 
laboratory  studies. 


3^2  Preliminary  Considerations 
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costs.  Proximity  to  the 
assist  the  field  operation; 
meteorological  station  will 
suitable  set  of  data  readily 
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gen  eral. 

while,  the 

presence  of 

a  nearby 

provide  the 

researcher 

wi  th  a 

available. 

3. 3  Location  and  Surf icial  Geology  of  the  Study  Area 


above 


field 

site  meeting 

most  of  the  criteria 

mentioned 

was 

selected 

near  Taber, 

Alberta 

Sec. 1 3- 

Twp.10-  Rg. 17- 

W4M)  ,  approximately 

50  km  east 

Canada 

Agriculture 

Lethbridge  Research 

S  ta  tion , 

ge,  Alberta  (Fig.  3- 

1).  The  study  area  is 

part  of  a 

large  terrace  above  the  escarpment  of  the  Oldman  River  (Fig. 
3-1)  It  is  entirely  underlain  by  glacial  till.  This  unit 
ranges  in  thickness  from  several  centimeters  to  more  than 
nine  meters.  In  the  immediate  vicinity  of  the  study  area, 
this  till  unit  is  found  at  a  depth  of  approximately  5.5 
meters.  Two  types  of  material  overlie  the  till;  (a)  alluvial 
sand  and  gravel  and  (b)  aeolian  sand.  In  the  study  area,  the 
deposits  overlying  the  till  consist  of  sandy  material  that 
is  classified  texturally  as  sandy  loam  to  sand.  Some 
discontinuous  thin  silt  or  clay  loam  lenses  are  interbedded 
with  the  sand. 


The  water  table  depth  is  variable,  ranging  from  2  to  3 
meters  below  the  surface  of  the  study  area.  In  general,  the 
lack  of  potholes  indicates  good  internal  drainage. 


* 
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Figure  3-1.  Location  of  study  area 


3.4  The  Field  Sites 


4  7 


Two  sites,  which  were  located  approximately  150  meters 
apart,  were  selected  for  detailed  study.  An  area 
approximately  9x9  meters  square  around  each  site  was  fenced 
to  protect  them  from  cattle  grazing  nearby. 

None  of  the  two  sites  was,  however,  entirely 
satisfactory.  The  extreme  sandiness  of  one  of  the  sites 
(Site  1)  severely  restricted  the  range  of  field-measured 
hydraulic  head  distribution  and  limited  the  eff ec ti vene ss  of 
laboratory  support  data  on  the  range  of  water  characteristic 
and  hydraulic  conductivity  functions.  The  second  site  (Site 
2)  ,  although  having  the  more  favorable  sandy  loam  textures, 
was  characterized  by  a  relatively  high  degree  of 
heterogeneity  of  the  subsurface  profile.  In  analysing  the 
field  data,  however,  more  emphasis  was  placed  on  data  from 
the  second  site. 

3._5  Weather  Observations 

The  Taber  area  is  generally  characterized  by  hot  and 
dry  summers  and  cold  winters  (Fig.  3-2).  A  small  weather 
station  is  located  approximately  five  kilometers  from  the 


field  sites 
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Meteorology 


Isohyet,  mean  annual  precipitation  in  cm 

Meteorological  station  . 

Standard  rain  gauge  only  . 


a 

a 


Precipitation  data: 

Mean  annual  precipitation  in  cm  . 

Year  of  commencement  of  observations  . 

Number  of  years  of  record  averaged  to  determine  the  mean  annual  precipitation  figure 


43.8 


1936 


10-24 


Mean  monthly  potential  evapotranspiration* . 

Mean  monthly  precipitation  . 

Period  when  surface  is  usually  snow  covered . 

Period  with  mean  daily  temperature  below  freezeing  (0°C) . 

Figure  indicates  percentage  of  mean  annual  precipitation  falling  as  rain 


Sources  of  data: 

Climatic  Maps  of  Alberta  (Longley,  1968);  Monthly  Record 
of  Meteorological  Observations  (Meteorological  Branch,  Canada 
Department  of  Transport);  and  Provisional  Evaporation 
Maps  of  Canada,  Department  of  Transport  Circular  4531 
(Bruce  and  Weisman,  1967). 


’Evapotranspiration  estimates  are  based  on  the  Provisional  Evaporation  Maps  of  Canada. 


(from:  Tokarsky,  1974) 


Figure  3-2.  Meteorology  of  study  area 
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3«_6  Measurement  of  Basic  Parameters  and  Equipment  Required 
for  the  Present  Study; 

In  the  following  paragraphs,  a  brief  description  of  the 
parameters  required  for  coupled  non-isothermal  field  studies 
will  be  presented  ;  for  instance,  types  of  materials, 
hydraulic  potential,  water  content,  water  characteristic 
functions,  hydraulic  conductivity,  density  and  porosity, 
temperature,  thermal  conductivity,  heat  capacity  and  water 
and  heat  fluxes.  Also  the  equipment  and  methods  of 
installation  used  to  measure  the  above  mentioned  parameters 
will  be  briefly  described. 

Hydrometeorological  equipment  at  the  study  area 
consisted  of  a  hygr other mograph  placed  in  an  instrument 
shelter,  a  Gen  atmometer  and  a  standard  rain  gauge.  This 
instrumentation  was  installed  at  Site  1.  Each  of  the  two 
sites  was  equipped  with  the  instrumentation  depicted 
diagrammat ica lly  in  Fig.  3-3.  Three  piezometers  were 
installed  at  Site  1  (la,  1b  and  1c)  and  two  at  Site  2  (2a 
and  2b)  in  June  1975.  The  rest  of  the  instrumentation  was 
installed  during  June  1976.  All  lead  wires  protruding  from 


the  ground 

surface 

were  protected 

by  passing  them  through 

rubber  rubes. 

the  ou 

ter  ends  of 

which  were 

covered  by 

plastic  bags 

.  For 

the  collection 

of  field  data  daily  or 

weekly  visits 

to  the 

sites  were  made 

at  irregular 

inte  rvals 

during  the  summer-early  fall  1976  measurement  period. 


hereafter  referred  to  as  the  1976  measurement  period 
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6  neutron  access  tube 

7  tensiometers 

8  hydrometeorological  equipment 


Figure  3-3.  Site  instrumentation  diagram 
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3.6.1  Soil  Ty  pes 


Both  sites  were  sampled  during  the  summers  of  1975  and 
1976  using  a  Stirling  auger  drill  and  a  hollow-stem  auger. 
Depth  intervals  of  approximately  30  cm  --  with  a  finer 
spacing  for  the  upper  layers  of  the  profile  —  were  sampled 
down  to  485  cm  depth.  Soil  texture  was  determined  in  the 
laboratory  using  the  hydrometer  method  (Day,  1965). 

3.6.2  Hydraulic  Potential 

3 . 6 . 2. 1  Piezo  meters 


Piezometers  we 

re  installed 

in 

order 

to  me 

asur  e 

saturated  hydraulic 

potentials  at 

the 

si tes. 

They 

were 

constructed  from  3. 

8  cm  I.  D.  PVC 

pipe. 

The  lower  30 

cm  of 

each  pipe  were  hacksaw  slotted  and  wrapped  with  fiberglass 
tape  to  form  an  intake,  and  the  lower  end  of  the  pipe  was 
capped.  The  piezometers  were  emplaced  by  lowering  them 
inside  of  the  15  cm  hollow-stem  auger  after  the  center-bit 
was  removed.  A  mixture  of  sand  from  the  hole  and  pea  gravel 
was  poured  around  the  outside  of  the  PVC  pipe  to  cover  the 
slotted  interval.  The  hole  was  then  filled  with  the 
excavated  material  and  the  surface  around  the  protruding  end 
of  the  pipe  was  sealed  with  bentonite.  The  upper  end  of  the 
piezometer  was  capped  using  a  PVC  fitting  with  a  small  hole 
drilled  in  it  for  pressure  equalization. 

Two  piezometers 


(Id  and  2d)  were  neither  slotted  nor 
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Piezometers  were  installed  at  the  following  approximate 
depths:  565  (1b),  435  (Id),  405  (1c)  and  350  (la)  cm  at  Site 

1  and  475  (2b),  440  (2d)  and  390  (2a)  cm  at  Site  2.  After 

installation,  the  piezometers  were  developed  by  repeated 
bailing. 


3. 6.2.2  Tensiometers 


Although  the 
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pressure  head  is 
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All  tensiometers  were  tested  in  the  laboratory  prior  to 
installation  in  the  field  and  the  porous  cup  conductance  and 
response  time  constant  for  most  of  them  were  determined 
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according  to  the  procedures  ou 
The  average  response  time  cons 
minutes.  Scales  directly  reading 
each  mercury  manometer  tensiomete 
reading. 

For  field  installation,  a  ho 
the  appropriate  depth  using  an 
King  tube.  Before  inserting  the  t 
loose  soil  from  the  hole  and 
poured  into  the  hole  to  insure  wa 
tensiometer  cup  and  the  soil.  T 
into  the  hole  with  a  twisting  dow 
surface  around  the  tensiometer 
amount  of  bentonite.  Tensiomet 
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icate  at  both  sites  to  dept 

90, 

120  and  180  cm 

in 

a  semicircl 
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{S ect ion  3.6. 

3)  . 

After  insta 

in 

each  nest  were 

filled  with  de 
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mercury  pot  was  filled. 

In  order  to  minim 

ize  the  inf 
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operation  of 
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field  instruments 
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taken  as 

possible  before 
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sun  had  w 

temperature  differ 
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from  that  o 

tlined  in  Danielson  (1  956}  . 
tant  was  approximately  2 
in  cm  H20  were  attached  to 
r  and  adjusted  for  the  zero 

le  was  made  in  the  soil  to 
approximately  2.5  cm  I.  D. 
ensiometer,  a  handful  of 
a  small  amount  of  water  was 
ter  continuity  betwen  the 
he  tensiometer  was  inserted 
nward  motion  and  the  soil 
was  sealed  with  a  small 
ers  were  installed  in 
hs  of  approximately  30,  60, 
e  around  the  neutron  access 
llation,  all  tensiometers 
aerated  distilled  water  and 

luence  of  temperature  on 
s,  the  first  reading  of  the 
early  in  the  morning  as 
armed  the  tensiometers  to  a 
f  the  soil. 
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3, 6. 2. 3  Soil  Psychrome ters 

Because  the  upper  limit  of  the  tensiometer  range  i 


about 

-0.8 

at  mo 

spheres, 
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potential  at 
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v| J  =  (RT/Mg)  In  h  (3-1 

Because  the  vapor  pressure  depends  on  both  the  matric  (o 
negative  pressure  head)  and  osmotic  components  of  the  soi 
water  potential,  psychrometers  measure  both  components  whil 
tensiometers  only  measure  the  matric  component  of  soil  wate 
potential. 

The  psychrometers  used  for  this  study  were  of  th 
Peltier  type.  They  were  constructed  following  procedure 
outlined  in  Kiebe  et  al.  (1971)  and  Lopushinsky  and  Kloc 
(1970),  except  that  ready-made  Chromel-Co nstantan  (0.025  m 
in  diameter  from  Omega  Engineering  Inc.)  thermocoupl 
junctions  were  used.  The  thermocouple  junction  was  protects 
by  a  cup-shaped  ceramic  bulb,  constructed  using  method 
described  by  Lopushinsky  and  Klock  (1970).  These  ce rami 
bulbs  were  checked  for  their  air-entry  value  by  means  of 
hanging  water  column.  The  ceramic  cup  was  considers 
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satisfactory  if  the  air-entry  value  was  in  the  range  of  at 
least  200  cm  H^O.  Prior  to  the  insertion  of  the  porous  cup, 
thermocouple  junctions  were  cleaned  carefully  using  an 
ultrasonic  device.  All  soil  psychrometer s  were  calibrated 
under  controlled  temperature  conditions  with  KCL  salt 
solutions  for  which  the  vapor  pressure  was  known.  Fig.  3-4 
shows  two  of  the  calibration  curves  obtained.  A  number  of 
ready-made  Wescor  Inc.  psychromet ers  were  also  used. 
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3.6.3  Water  Table  and  Water  Content 


Water  table  wells  were  constructed 
throughout  most  of  its  length,  wrapping 


by  slotting  PVC  pipe 
it  with  fiberglass 


and  emplacing  it  in  boreholes.  The  volumetric  water  content 
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Psychrometer  75-20 


-5  -10  -20 


Water  potential  - bars 


Figure  3-4.  Psychrometer  calibration  curves 


from  a  depth  of  approximately  20  cm  below  the  ground  su 
down  to  the  water  table  was  measured  using  the  ne 
scattering  method  (Gardner,  1965;  Long  and  French,  1967 
each  site,  an  aluminum  access  tube,  approximately  4  m 
long,  4.0  cm  I.  D.  and  sealed  at  its  bottom,  was  insert 
an  augered  hole.  A  Nuclear-Chicag c  model  5810  neutron 
probe  in  conjunction  with  a  portable  scaler  model  592 
used  to  take  three  one-minute  counts  for  each  15  or  3 
depth  interval.  Standard  counts  were  taken  before  and 
water  content  measurements.  The  ratio  of  the  average 
to  the  average  standard  count  was  used  in  a  1 
regression  equation  provided  by  the  manufacturer  to  est 
the  water  content  on  a  volumetric  basis. 

Prior  to  field  use,  the  neutron  probe  was  calib 
using  a  sand  with  a  composition  similar  to  that  at  the 
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3.6.4  Water  Characteristic  Curves 
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Fiaure  3-5.  Neutron  calibration  curve 


59 


(water  characteristic  curv 
simultaneous  field  measurem 
the  tensiometers.  However, 
contents  and  pressure  head 
period,  characteristic  cu 
laboratory  in  order  to  chec 


the 

ran 

ge 

o 

f 

me  as 

ureme 

n  t  s 

• 

Fo 

ur 

se 

para 

te  m 

ethod 

soil 

sa 

mple 

s 

in  t 

he  la 

borat 

char 

act 

er 

is 

ti 

c  c 

ur  ves 

o 

ve 

tens 

ion 

» 

F 

or 

ten 

sions 

ab 

ov 

techniq 

ue 

s 

i 

nvol 

ving 

sa  tur 

(Dan 

iel 

so 

n. 

1 

956) 

.  For 

me 

as 

te  ns 

ion 

ra 

ng 

e. 

a  pr 

oce 

du 

1965 

)  i 

nv 

ol 

vi 

ng  a 

pres 

sur 

e 

the 

t  e 

ns 

io 

n 

range  of 

100 

t 

cera 

mic 

P 

or 

ou 

s  pi 

a  tes 

wer 

e 

A 

sat 

ur 

at 

io 

n-ca 

pilla 

ry 

experim 

en 

t 

(L 

alib 

erte 

et. 

al 

the 

wat 

er 

c 

ha 

ract 

erist 

ic 

fu 

mb. 

Da 

ta 

fr 

om 

gravi 

met 

ri 

convert 

ed 

t 

o 

volu 

metric  m 

oi 

the 

bul 

k 

de 

ns 

ity 

of  th 

e  s 

oi 

each 

site  was  obtained 

from 

sing 

the  neutron  probe 

and 

the 

limited  range  of 

water 

rved 

during  the  measurement 

were 

determined  in 

the 

field 

results  and  to  e 

xtend 

e 

u 

tili 

zed 

fo 

r 

test 

in 

g 

the 

n 

order 

to 

ob 

tain 

wa 

ter 

su 

ff 

icie 

nt  ly 

wi 

de  r 

an 

ge 

of 

bars 

,  v 

apor 

ea 

-i. 

uili 

br 

at 

ion 

2S 

°+ 

solu 

tion 

s 

we 

re  e 

mp 

loyed 

ts 

in 

the 

1 

to 

15 

bar 

sc 

ri 

bed 

by 

Rich 

ards 

( 

IS 

54, 

ne 

a 

ppar 

at  us 

w 

as 

use 

d. 

For 

0 

mb 

/ 

essu 

re 

c 

ooke 

rs 

w 

it  h 

ed 

(Rich 

ards 

/ 

19 

54, 

1 

96 

5)  . 

re 

( 

or 

pr 

es 

sure 

he 

ad) 

6) 

w 

as  e 

mplo 

yed 

to 

me 

as 

ur  e 

for 

ten 

sion 

s 

le 

ss  t 

ha 

n 

100 

St 

ur 

e  d 

eter 

mi 

na 

tion 

s 

w 

ere 

cont 

en  ts 

by 

mu 

It 

iply 

in 

g 

by 

les. 


e)  at 

en ts  u 

due  to 

s  obse 

r  ves 

k  the 

s  wer 

or  y  i 

r  a 

e  15 

ated  H 

;uremen 

re  de 

membra 

o  1,00 

utiliz 

pressu 

.  ,  196 

notion 

c  moi 

s  tur  e 

1  samp 


60 


3.6.5  Hydraulic  Conductivity 


3. 6. 5. 1  Saturated  Hydraulic  Conductivity 


The  saturated  hydraulic  conductivity  was  measur 
the  field  using  the  piezometer  method  (Kirkham, 
Hvorslev,  1951).  When  the  inlet  section  of  a  piezomet 
constructed  by  slotting  the  casing,  hydraulic  conduct 
is  measured  essentially  in  a  horizontal  direction  while 
only  the  end  of  the  casing  available  as  an  intake,  hudr 
conductivity  is  determined  in  a  predominantly  ver 
direction. 

The  hydraulic  conductivity  is  determined  by  the  us 
the  following  formula  based  on  Darcy's  equation  (Kir 
1945) 

K  =  [*R2  In  (y,  /y4  )  ]/[  S  (t4-t<  )  ] 
or,  using  Hvorslev's  (1951)  definition  of  basic  time  la 

K  =  /IR2/(ST') 

where  R  is  the  radius  of  the  piezometer; 

y1  is  the  distance  from  the  water  table  to  the 
level  in  the  piezometer  at  time  ti  ; 
y  is  the  distance  from  the  water  table  to  the 
level  in  the  piezometer  at  time  t2; 

S  is  a  shape  factor  with  dimensions  of  length 

depends  on  the  geometry  of  the  inta 
the  piezometer. 
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Briefly,  the  determination  procedure  was  as  follows.  An 
average  characteristic  (desaturation)  curve  (  v|/ vs  Sw  or  9  ) 
obtained  from  three  soil  samples  was  derived  and  the 
relationship  between  hydraulic  conductivity  and  pressure 
head  was  determined.  The  methods  and  equipment  employed  were 
those  described  in  Brooks  and  Corey  (1964) ,  Anat  et  al. 
(1965),  Laliberte  et  al.  (1  966),  van  Schaik  and  Laliberte 
(1968)  and  van  Schaik  and  Graham  (1969).  From  the 
desaturation  curve,  the  residual  saturation  Sr  was 
determined  by  a  trial  and  error  procedure  (Brooks  and  Corey, 
1964).  The  effective  saturation  Sa  ,  defined  by  Brooks  and 
Corey  (  1  964)  as 

Se  =  (Sv/-Sr)/(1-Sr)  (3-6) 

is  linear  for  Se  less  than  approximately  0.8  when  plotted 
versus  on  a  log-log  paper.  The  slope  of  this  straight  line 
portion  is  called  the  pore-size  distribution  index 
(Brooks  and  Corey,  1964).  The  intercept  of  the  straight  line 
at  Se=1.0  is  called  the  bubbling  or  air  entry  pressure  vk, 
which  corresponds  very  closely  to  the  minimum  pressure  head 
at  which  air  begins  to  enter  a  saturated  porous  medium  and 
the  gas  phase  first  becomes  continuous.  The  relative 
hydraulic  conductivity  Kr  (=K/Ksa^  )  may  then  be  calculated 
as 

Ky_  =  Se  <2+3^  >/^7  (3-7) 

or 


Ky  =  ('Vfc/'p)n 


where  2+3iV. 


for 


(3-8) 
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Accordingly,  by  combining  the  experimental  desaturation 
data  with  the  calculated  value  of  residual  saturation,  the 
relationship  between  ^  and  $  from  saturation  to  residual 
saturation  is  obtained.  When  the  experimentally  determined 
Ky.  vs  relationship  is  plotted  on  a  log-log  paper,  a 

straight  line  results  for  ,  as  is  implied  by  (3-8) 

above.  This  line  has  a  slope  of  — and  intercepts  the  ^-axis 
at  the  air-entry  pressure  v|^,  where  Kr  =1.0.  This  straight 
line  may  be  extrapolated  to  higher  absolute  values  of 
Thus,  by  combining  this  Ky.  vs  relation  with  the  vs  $ 
relationship  for  the  same  soil,  the  hydraulic  conductivity 
as  a  function  of  water  content  from  residual  saturation  to 
saturation  is  obtained. 
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ated  water  flux  qw/^>  is  known 
the  isothermal  and  non-iso thermal 
(2-38)  in  one-dimensional  form  may  be 
Thus,  using  the  non-isother raal  water 
8)  in  combination  with  (2-37),  the 
conductivity  may  be  calculated  as 
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follows 
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3.6.6  Density  and  Porosity 


Bulk  density  was  calculated  from  cores  of  k 
obtained  in  the  field  to  depths  of  approximate 
Alternatively,  a  Cs-gamma  density  probe  (Mode 
also  used  to  measure  bulk  densities  at  various  de 

Particle  densities  were  determined  using  the 
method  (Blake,  1965).  Once  bulk  and  particle 
known,  porosity  can  be  determined  from  th 

relationship: 

t  =  1 

where  ^  is  bulk  density  and 
^  is  particle  density. 
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Air  temperatures  were  measured  using  a  hygro thermo  graph 
which  provided  continuous  temperature  and  air  humidity 
recordings  over  weekly  periods.  The  hygr ot hermograph ,  which 
was  also  checked  for  accuracy  in  the  laboratory,  was  placed 
in  an  instrument  shelter  approximately  1.3  meters  above  the 
ground  surface. 

3^6^ 8  Thermal  Conductivity 
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The  thermal  conductivity  probes  used  for  this  study 
were  similar  to  those  described  by  Slusarchuk  and  Foulger 
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(1  973).  1  The  following  modi 
thermistor  (GB35L  1-Fenwal  E 
was  potted  in  a  hole  in 
tube.  Teflon-coated  constan 
plus  3  mil  teflon  wall) 
manner  around  the  plastic  t 
heating  element  approxima 
wires  (  2  thermistor  leads 
connected  to  a  four-wire 
ring,  3.8  cm  in  diameter,  w 
and  small  protruding  rod 
attached  near  the  top  of  th 
tube.  A  plastic  cap  was 
ring  to  seal  and  protec 
junctions.  An  arch-shaped 
attached  to  the  ring.  A  Ion 
arch  for  recovery  purposes 
checked  for  leaks  by  ins 
long  pipe  filled  with  water 


fications  were  employed.  0 
lectronics,  "Thermistor  M 
the  wall  of  a  4.85  mm  OD 
tan  wire  (  of  0.127  mm  d 
was  wound  tightly  in  a 
ube  using  a  lathe  to 
tely  33  cm  long.  The  fo 
and  2  constantan  leads 
coaxial  cable.  A  stainles 
ith  a  groove  near  its  pe 
s  from  two  opposite  sid 
e  8.0  mm  O.D.  stainless 
then  glued  into  the  groove 
t  the  lead  wire-coaxial 
thin  stainless  steel 
g  steel  wire  was  connected 
(Fig.  3-6).  The  probe 
erting  them  for  several  da 
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In  order  to  insert  and  recover  the  probes  at  depths 
below  the  water  table,  a  device  was  designed  (Fig.  3-7)  to 
fit  and  lock  on  the  two  protruding  rods  of  the  ring 
mentioned  above.  This  device  was  connected  to  steel  rods 


details  on  the  operation  and  construction  of  thermal 
conductivity  prooes  may  also  be  found  in  Wechsler  (1  966)  ; 
Janse  and  Borel  (1965);  Jackson  and  Taylor  (1965);  Beck  et 
al.  (1971),  Fritton  et  al.  (1974),  among  others. 
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Figure  3-6.  Thermal  conductivity  probe 
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which  fit  together  to  facilitate  installation  in  bore 
of  any  desired  depth.  Aluminum  conduit  tubes,  approxim 
5  cm  in  diameter,  were  inserted  in  drilled  holes  at  d 
of  approximately  305  and  340  cm  in  Site  1,  and  375  c 
440  cm  in  Site  2.  The  conductivity  probes  were  emplaced 
the  soil  at  the  bottom  of  the  conduit  tubes  by  using 
device  described  above.  At  shallower  depths  of  approxim 
6,  30  and  60  cm  the  probes  were  simply  inserted  horizon 
into  the  walls  of  the  trench  dug  at  each  site.  In •all  c 
only  the  heating  curve  was  considered.  The  duration  of 
experiment  in  the  field  was  approximately  12  minutes. 
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Figure  3-7.  Thermal  conductivity  probe  insertion 
and  recovery  apparatus 


71 

plotted  curve.  Comparison  of  the  mean  value  of  thermal 
conductivity  for  the  c oarse-grain ed  dry  Ottawa  sand  obtained 
from  five  probes  ( ^  =7. 21 4X 1 0~4  cal/cm/sec/C)  by  Slusarchuk 
and  Foulger  (1973)  as  well  as  of  their  guarded  hot  plate 
value  (7.  262X  10-*  cal/cm/sec/C)  --  which  was  treated  as 
their  thermal  conductivity  standard  --  with  the  mean  value 
{  ( 7 .  1 1 4  ±0  .  1 68  )  XI  0-<v  cal/cm/sec/C}  obtained  from  our  ten 
probes  under  the  same  conditions,  showed  a  satisfactory 
correspondence  (Table  3-1).  For  calibration  purposes,  the 
ratios  of  our  air-dry  and  saturated  Ottawa  sand  thermal 
conductivities  for  each  probe  to  those  of  Slusarchuk  and 
Foulger  (1973)  using  the  guarded  hot  plate  apparatus  were 
joined  by  straight  lines,  assuming  a  linear  relationship 
between  thermal  conductivity  and  water  content.  Therefore, 
at  intermediate  water  contents  between  air  dryness  and 
saturation,  the  calculated  thermal  conductivity  may  be 
divided  by  the  above  mentioned  ratio  at  the  appropriate 
water  content  to  give  a  standardized  thermal  conductivity 
value. 


Several  empirical  (Kersten,  1949)  and  semi-empirical 
equations  (de  Vries,  1963)  have  been  developed  for 
calculating^.  However,  as  van  Booyen  and  Winterkorn  (  1957) 
indicated,  none  of  those  equations  was  able  to  give 
consistently  dependable  results.  The  thermal  conductivity 
probe  method  avoids  all  problems  involved  with  a  theoretical 
calculation  of  thermal  conductivity  (Kimball  et  al. ,  1976, 


Hadas,  1977) 


&  **r  c,  ■  tii  ifri  "j 
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Figure  3-8.  Field  equipment  for  thermal  conductivity 


measurements 


time - ►sec 
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Measured  thermal  conductivity  of  air-dry  20/30  mesh  Ottawa  sand 
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3.6.9  Heat  ca pac ity 

The  volumetric  heat  capacity  of  a  soil  (de  Vries,  1963) 
can  be  found  by  adding  the  heat  capacities  of  the  mineral 
matter  (C^)  t  the  organic  matter  (C0),  the  water  (Cw)  and 
the  air  constituents  in  a  unit  volume  of  porous  medium.  The 
heat  capacity  of  air  is  very  small  and  can  be  neglected. 
Thus,  the  volumetric  heat  capacity  of  a  soil  is  given  as 

c  =  +9cw  (3-14) 
where  C  is  the  volumetric  heat  capacity  (cal/cm3/C)  of  a 
moist  porous  medium  and  x*,*, ,  x0  are  the  volumetric  fractions 
of  mineral  matter  and  organic  matter  respectively,  and  0  is 
the  volumetric  water  content.  de  Vries  (1963)  gives  the 
following  average  values  of  heat  capacities:  CW)n=0.46 
cal/cm3/C  and  C0  =0.60  cal/cm3/C  .  The  volumetric  heat 
capacity  of  water  is  unity.  Therefore,  the  above  equation 
becomes 
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In  order  to  obtain  the  volumetric  percentages  of 
organic  and  mineral  fractions  of  the  various  soil  layers,  an 
average  organic  matter  density  of  1.3  gm/cm3  (de  Vries, 
1S63)  was  assumed.  The  average  mineral  particle  density  was 
found  to  be  2.64  gm/cm3  in  the  top  60  cm  of  the  profile. 
This  value  was  assigned  to  the  entire  subsurface  profile. 
Using  these  density  values  together  with  the  field-measured 
bulk,  densities,  the  organic  and  mineral  matter  fractions  x0 


and  were  calculated  as  follows: 


(3-16) 

(3-17) 

mea  sured 


Xo  =  <*'  (  ) 

XWW,  =  (1-<**)  ?b/^5 

where  o(‘  is  the  organic  matter  fraction  by  weight, 
experimentally.  Values  of  x0  and  at  each  depth  in  the 

profile  remain  constant  for  non-def ormable  media.  However, 
because  the  volumetric  water  percentage  9  varies  in  time, 
the  volumetric  heat  capacity  (3-14)  varies  as  a  function  of 
time. 


3. 6. 10  Heat  and  Hater  Fluxes 


3. 6.  10.  1  Heat  Flux 


The  heat  fluxes  at  different  depths  in  the  subsurface 
profile  were  computed  by  the  temperature  integral  or 
calorimetric  method  (Portman,  1954;  Lettau  and  Davidson, 
1957;  Tanner,  1963;  Kimball  and  Jackson,  1975).  In  this 
method,  the  heat  flux  is  calculated  from  the  change  of  the 
heat  content  or  storage  in  the  subsurface  profile  over  a 
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given  time  interval-  The  equation  for  the  heat 
computation  can  be  derived  from  the  heat  balance  equ 
(2-62) ,  which  in  one-dimensional  form  is 

C  5T/9t  =  ( 

Multiplying  both  sides  of  (3-18)  by  dz  and  integrating 
depth  z  to  depth  z+dz  results  in 

=  gyjxvd*  +  Jc  3T/at  dz  ( 
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only  determine  heat  flow  by  conduction  because  the  met 
impervious  to  liquid  and  gas.  Therefore,  heat  transfe 
water  liquid  and  vapor  movement  cannot  be  determined 
same  limitation  is  apparent  with  the  null-alignment  m 
(Kimball  and  Jackson,  1975),  where  null  points  1 
temperature  gradient  are  used  to  provide  zero  heat  flux 
known  depths  in  the  profile. 


er  is 
r  by 
.  The 
ethod 
n  the 
es  a  t 


3.6.10.2  Water  Flux 


The  total  water  flux 

in 

the 

unsat ur at 

ed 

zone  can 

be 

measured  in  the 

same 

way 

as 

the  heat 

flux 

;  that  is. 

by 

accounting  for  the 

change 

of 
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he  subsurf 

ace 

profile  over  a  given  time  interval.  If  a  one-dimensional 


water 
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whole  profile 
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calculated  u 
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the 

above 

equation  (Richa 

rds 

e t  al.  , 

19 

56; 

Rose  et  al. ,  1965;  Rose  and  Krishnan,  1967;  van  Ravel  et 
al.,  1  968;  Flocker  et  al.  ,  1  968;  Hillel  et  al.  ,  1972  ; 


Nielsen 

et  al. , 

1  973) 

.  However,  a  reference 

water 

flux 

at  a 

certain 

dept  h 

may 

not  be  easy  to  determine. 

The 

ma  jo  r 

problem 

with 

this 

approach  is  that 

the 

hy  dr 

aulic 

conductivity  function  must  be  determined.  The  procedure  I 
have  adopted  is  the  following:  (1)  sample  soil  from  any 
layer  of  the  subsurface  profile  and  determine  the  hydraulic 
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co nduct iv ity- wate r  content  function  in  the  laboratory;  (2) 
use  the  non-isother mal  flux  equation  (2-38)  in 
one-dimensional  form  to  calculate  the  water  flux  at  the 
bottom  of  this  layer,  the  K  function  being  known.  Other 
relevant  parameters,  such  as  the  depth  distribution  of 
hydraulic  head,  are  already  known. 


The  procedure  followed  by  Nielsen  et  al. ,  1964,  1973; 
Davidson  et  al. ,  1969,  among  others  —  which  involves 
covering  the  ground  surface  with  plastic  sheets  and  assuming 
zero  water  flux  (no  evaporation)  at  the  surface  --  may  not 
be  appropriate  under  conditions  not  near  saturation. 
Suppression  of  evaporation  will  increase  the  subsurface 
temperature  regime  and  change  the  water  vapor  flux 
contribution  to  total  water  flow.  The  method  used  by  Rose 
(1968c)  for  calculating  bare  soil  evaporation  by  application 
of  the  principle  of  water  conservation  (3-9)  in  a  volume  of 
surface  soil,  in  combination  with  the  non-isothermal  water 
flux  equation  (2-38)  in  one-d imensiona 1  form,  may  be 
preferable.  Application  at  the  field  sites  of  the  method 
proposed  by  Brutsaert  (1975)  for  calculating  evaporative 
fluxes  based  on  non-dimensionless  parameters  (similarity 
concepts)  resulted  in  unacceptably  high  evaporative  fluxes. 
Therefore,  this  approach  was  not  taken  further.  With  regard 
to  the  saturated  water  flux,  this  may  be  calculated  from  the 
heat  flux  equation  as  suggested  in  (2-72).  (An  evaluation  of 
such  a  method  is  presented  in  section  3.8). 


. 
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The  potential  evaporative  flux  from  a  saturated  surface 
was  measured  at  the  field  sites  using  a  Gen  atmometer 
(Mukammal,  1961).  This  atmometer  consists  of  a  flat  black 
porous  plate  (Bellani  plate),  7.5  cm  in  diameter,  kept  moist 
by  distilled  water  held  in  a  tube  which  is  connected  to  a 
supply  reservoir. 

3.7  Experimental  Results  and  Discussion 


3.7.1  Soil  Ty_pes  and  Analyses 


Site 

and  3-11)  o 
Underly ing 
below  appro 
lens  was 
variation  w 
Site  1  sand 
profile  i 
cond  uct ivit 
probably  d 
porous  mate 

Site  2 
Sit  e  1 .  T 
laterally  a 
(Fig.  3-11 
soil  textur 


1 

is 

com 

prised 

of  unif 

or 

m 

me 

di 

urn 

sand 

(Figs. 

3-10 

■f 

an 

ave 

rage 

P 

article 

de 

ns 

it 

y 

of 

2.  65 

42  gm/cm3. 

t 

he 

sa 

nd 

is 

a  cla 

y 

lo 

am 

t 

ill  foun 

d  at  depths 

xi 

mately 

5.5 

m 

.  Howeve 

r. 

a 

n 

i 

so 

lated 

clay 

loam 

sa 

mpled 

at 

a 

430-46 

0 

c 

m 

d 

ep 

th  in 

ter val 

.  The 

it 

h 

depth  of 

o 

ther  ph 

ys 

ic 

al 

pr 

oper ti 

es  of 

the 

i 

s 

indi 

cate 

d 

in  Fig. 

3- 

10 

• 

Al 

though  t 

he  tex 

tural 

s 

near 

iy 

uniform , 

th 

e 

s 

at 

urated 

hy  dr 

aulic 

Y 

is 

var 

iabl 

e 

(Fig. 

3- 

10 

). 

Th 

is  si 

t uatio 

n  is 

ue 

to  d 

if  f  e 

re 

nces  in 

th 

e 

CO 

mpac 

tion  s 

tate  o 

f  the 

ri 

al 

,  as 

well 

as  to  ex 

pe 

ri 

me 

nt 

al 

error 

s. 

h 

as 

a  more 

CO 

mplicate 

d 

ge 

olog 

ic 

al  se 

tting 

than 

he 

subs 

urf  a 

ce 

units  a 

re 

t 

ex 

tu 

ra 

lly  va 

riable 

bot  a 

nd 

vert 

ically 

,  w  it  h 

s 

an 

dv 

lo 

am  pr 

edomin 

ating 

). 

Fig. 

3- 

12 

shows 

th 

e 

va 

ri 

at 

ion  wi 

th  d  ep 

th  of 

es 

r 

as  w 

ell 

as 

other 

phys 

ic 

al 

proper 

ties 

of  a 

,i  <a-r;  (c 


►  gm/cc  K  - ►  cm /hr  %0M 


81 


0  50  100  150  200  250  300  350  400  450 


z - ►  cm 


Figure  3-10.  Physical  properties  of  Site  1 
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Figure  3-11.  Soil  texture  triangle  for  the  study  area 


sampled  borehole  at  Site  2.  A  correlation  between  perce 
of  sand  and  saturated  hydraulic  conductivity  is  app 
(Fig.  3-12).  The  average  particle  density  of  the  top  3 
of  the  soil  was  2.6281  gm/cm3. 


3.7.2  Hydraulic  Head 
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The  time  distribution  of  hydraulic  head  as  indicat 
tensiometers  and  piezometers  at  Site  1,  as  well  as  the 
table  fluctuations  during  a  nineteen-hour  measur 
period,  are  shown  in  Fig.  3-14.  The  advantage  of 
prolonged  measurements  is  that  the  diurnal  change 
several  physical  parameters  can  be  studied  in  deta 
assess  their  signif icance.  From  Fig.  3-14  it  can  be  de 
that  the  direction  of  flow  was  generally  downwards.  How 
the  direction  of  flow  was  not  always  constant  dependi 
whether  water  was  infiltrating  after  a  rainfall  or 
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Figure  3-12.  Physical  properties  of  Site  2 
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Figure  3-13.  Time  distribution  of  piezometric  water  levels 
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evaporating  during  an  extended  drying  period.  F 
the  period  of  measurements  depicted  in  Fig.  3- 
11-12,  1976)  followed  a  period  of  significant  rai 

3-13)  and  water  was  still  percolating  downwar 
toward  the  end  of  the  measurement  period,  the  hyd 
gradient  between  tensiometer  6  and  piezomete 
reversed,  indicating  the  completion  of  the 
process  and  a  return  to  evaporation-controlled  up 
The  piezometric  head  distribution  indicated  dow 
flow  at  shallow  depths  superimposed  on  a  deeper  u 
flow  condition  (Fig.  3-14).  This  saturated  flow  p 
consistent  throughout  the  1976  measurement  period 
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The  hydraulic  head  distribution  at  Site  2  was  more 
complicated  and  variable  than  at  Site  1.  The  significant 
degree  of  heterogeneity  of  the  subsurface  profile  of  Site  2 
contributed  to  this  hydraulic  head  pattern.  Fig.  3-15 
presents  a  twenty-four-hour  measurement  of  hydraulic  head, 
where  a  reversal  of  hydraulic  head  gradient  between 
tensiometers  4  and  6  is  observed  twice  during  that 
measurement  period.  However,  the  observed  hydraulic  head 
distribution  in  the  saturated  zone  with  an  upward  water  flow 
pattern  at  shallow  depths  superimposed  on  a  deeper  downward 
water  flow  pattern  was  consistent  throughout  the  1976 
measurement  period. 

Despite  the  care  taken  in  constructing  and  calibrating 
the  soil  psychrometers,  their  performance  in  the  field  was 


hydraulic  head 
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Figure  3-14.  Time  distribution  of  hydraulic  head  at  Site  1 
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not  satisf ac tor y ,  Normally,  soil  psychrometers  are  most 
useful  for  dry  soils  with  tensions  above  one  atmosphere 
because  of  the  higher  resolution  of  the  psychrome ter  output. 
However,  after  installation,  an  extended  wet  period  followed 
with  soil  water  tensions  usually  in  the  tensiometer  range. 
This  situation  reduced  the  effectiveness  of  psychrom  etric 
measurements.  Although  occasional  good  psychrometer  readings 
were  obtained,  they  were  not  consistent.  In  addition,  the 
calibration  of  the  psychrometers  left  in  the  field  probably 
changed  with  time  as  a  result  of  corrosion  of  the 
psychrometer  wire  (Wiebe  et  al. ,  1971).  It  was  considered 
impractical  to  remove  and  recalibrate  the  psychrometers 
during  the  course  of  the  field  measurement  period.  The 
procedure  of  encasing  thermocouple  psychr omet er s  in 
perforated  stainless  steel  tubes  (Moore  and  Caldwell,  1972) 
would  have  facilitated  removal,  recalibration  and 
reinsertion  in  the  field. 


3.7.3  Water  Contents 


A  typical  distribution  of  water  co 
obtained  by  using  the  neutron  prob 
3-16.  Due  to  the  very  uniform  grain  siz 
water  content  distrinution  indicated 
smooth  curve.  In  contrast,  the  distribu 
depth  at  Site  2  is  complicated  (Fig.  3- 
lower  hydraulic  conductivity  layers,  s 
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loam  (SL)  layers,  produces  the  observed  moisture  "bulges” 
just  above  those  low  conductivity  layers.  This  water  content 
distribution  with  depth  was  persistent  throughout  the  1976 
measurement  period. 
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Fiqure  3-16.  Time-depth  distribution  of  water  contents  at 
Site  1 
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Figure  3-17.  Time-deDth  distribution  of  water  contents  at 
Site  2 


(Figs.  3-18  and  A-1) 
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Figure  3-18.  Time-depth  distribution  of  temperature  at  Site  2 
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Figure  3-19.  Water  characteristic  curve  for  Site  1 
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Figure  3-20.  Water  characteristic  curve  for  Site  2  (0-45  cm) 


cm  H2O 


99 


_ 1 - 1 _ 

15  20  25 


%0 - * 


Figure  3-21.  Field  water  characteristic  curve  for  Site  2  (45-180  cm) 
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Figure  3-22.  Overall  characteristic  curve  for  Site  2 
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Figure  3-23.  Hysteresis  loop 
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tests  together  with  Ks<xt 

values  determined  in 

the 

labor  ator  y 

using  the  constant  head 

method.  In  general. 

the 

ag  re  emeu  t 

between  field  and 

laboratory  values 

of 

hy dr  aulic 

conductivity  and  their  correlation  with  soil 

tex 

ture  falls 

within  the  expected 

range  of  reliabi 

lity 

for  such 

measure  ment s. 


The  results  or  the  laboratory  hydraulic  conductivity 
measurements  are  presented  in  two  sets  of  curves  for  each 
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the 
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Al  th 

ough 

th 

e  a  bo 

ve 

me 

nt 

ioned 

method  of  hydraulic  conductivity  measurement  provides  a 
satisfactory  tool  for  extending  the  hydraulic  conductivity 
values  through  the  linear  portion  of  the  Kr (^)  relationship, 
caution  should  always  be  exercised  in  interpreting  and 
quantifying  such  extrapolated  results. 
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Figure  3-24.  «r  U)  relationship  for  Site  2  (0-20  cm) 
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Figure  3-25.  K  (e)  relationship  for  Site  2  (0-20  cm) 
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Figure  3-26.  Field-determined  K  (e)  curve  for  Site  2 
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3.7.7  Thermal  Conductivity  and  Heat  Capacity 

Figure  (3-27)  shows  a  typical  field  curve  obtained 
during  the  month  of  July  1976  from  Site  1  using  the  thermal 
conductivity  probe.  The  value  of  thermal  conductivity 
calculated  from  these  data  is  also  indicated.  The  points  on 
such  graphs  asymptotically  approach  a  straight  line 
somewhere  between  100  and  200  seconds  after  the  start  of  the 
test.  The  early  time  deviation  is  due  to  initial  lag  effects 
(Wechsler,  1966)  mainly  attributable  to  the  effects  of  probe 
diameter  and  conductivity,  and  contact  resistance  between 
the  test  medium  and  the  probe. 

Assuming  that  the  entire  subsurface  profile  is 
homogeneous  with  respect  to  thermal  conductivity,  a 
relationship  between  thermal  conductivity  and  water  concent 


M9)  is 

represented  in 

Fig.  3-28 

by  a  hand-drawn 

line 

passing 

thro  ugh 

the 

points.  The 

limited  range  of 

water 

contents 

over 

which 

the  thermal 

cond  uctivities 

were 

determined  in  the  field,  as  well  as  thermal  conductivity 
heterogeneities,  prevented  a  better  definition  of  the  A(0) 
relationshi  p. 

The  volumetric  heat  capacity  C (0 )  function,  calculated 
from  the  procedure  outlined  in  Section  3.6.9,  is  shown  in 
Fig.  3-29.  The  deviations  from  non-linearity  observed  in 
that  figure  are  due  to  the  variation  in  bulk  density  and 
organic  matter  content  (Fig.  3-12). 
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Figure  3-27.  Field-determined  thermal  conductivity 
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Figure  3-28.  X  (e)  relationship  for  Site  2 
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Figure  3-29.  C  (e)  relationship  for  Site  2 
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3. 8  Evaluation  of  Field  and  Laboratory  Studies 


Reliable  field  measurements  can  be  achieved  by  using 
properly  designed  and  calibrated  field  instrumentation  in 
combination  with  careful  field  installation  involving 
minimum  soil  disturbance  and  appropriate  spacing.  However, 
even  if  appropriate  precautions  are  taken,  the  field 
researcher  is  faced  with  many  difficulties  in  assuring  the 
long-term  accuracy  of  the  instrumentation.  For  example,  in 
the  case  of  the  mercury  tensiometers,  the  problems  of 
maintaining  the  associated  manometer  completely  full  of 
liquid  under  partial  vacuum,  time  lag  and  temperature 
influences  on  tensiometer  performance  cause  appreciable 
difficulties  in  obtaining  accurate  measurements  of  matric 
potentials  over  a  long  period  of  time.  Other  problems  with 
tensiometers  relate  to  difficulties  in  maintaining  constant 
physical  contact  with  the  soil,  especially  with  windy 
conditions;  disruptions  caused  by  repeated  servicing;  and 
plant  roots  which  tend  to  grow  around  the  base  of  the 
tensiometer. 

As  mentioned  earlier,  the  field  use  of  soil 
psycnro mete rs  was  not  satisfactory.  This  result  is  typical 
of  problems  reported  by  other  workers  (Maclean,  1974;  Mayo, 
1976).  Thus,  few  successful  field  measurements  with  soil 
psycnro me ters  are  reported  in  the  literature  (Enfield  and 
Hsieh,  1971;  Merril  and  Rawlins,  1 972  ;  Moore  and  Caldwell, 

major  shortcomings  of  this 


1972;  Wheeler  et  al. ,  1972) .  The 


, 


instrumentation 

are 

related 

to  the  delicacy 

of 

psycnrome tr ic 

measurements  them 

selves  and 

the 

f  re 

recalibration 

of  field 

psychrome 

ters  necessa 

ry  for 

acc 

measure  me  nt  s. 

The  variation  of 

water  co 

ntent  and 

temperat  ur 

usually  very 

small  under  normal 

conditions 

and  it 

is 

not  possible 

to  measure  these 

change  s. 

Only 

wh  en 

magnitude  of  the  changes  are  large,  such  as  in  the 
soil  zone  and  capillary  fringe,  can  these  change 
measured  with  reasonable  accuracy.  Difficulties  in  meas 
small  temperature  and  moisture  gradients  also  influenc 
accuracy  of  derived  data  which  are  functions  of 
variables.  Thus,  the  possibility  of  accurately  speci 
the  heat  flux  and  particularly  the  water  flux  at 
location  in  the  subsurface  profile  is  not  great 
addition,  the  small  magnitude  of  several  parameters  inv 
in  the  calculation  of  the  subsurface  flow  regime 
especially  of  the  thermal  diffusion  coefficients,  cou 
coefficients  and  conductivity  coefficients  add  to 
difficulties  in  applying  flow  theories. 


The  methodology  proposed  in 
make  the  field  endeavor  easier.  H 
of  relevant  quantities  often  po 
For  example,  I  found  that  calcula 
flux  from  the  heat  flux  equation 
Stallman  (1963)  and  Bredehoeft  a 
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questionable  for  field  situations. 


The  field  methodology  depends  strongly  on  the  accurate 
determination  of  water  and  especially  of  the  heat  fluxes  in 
the  subsurface  profile  (Taylor  and  Cary,  1965).  However, 
relatively  small  errors  in  the  calculation  of  heat  and  water 
fluxes  propagate  and  amplify  in  subsequent  transport 
coefficient  calculations.  For  example,  let  HF  be  the  ratio 
of  the  heat  flux  calculated  using  the  temperature  integral 
method  (Section  3.6.10.1)  over  the  heat  flux  calculated 
using  (2-58)  or  (2-60).  Also  let  DC  be  the  ratio  of  the 
water  vapor  diffusion  coefficient  Dqv  --  (2-71)  neglecting 
the  convective  term  --  estimated  using  the  heat  flux 
calculated  by  the  temperature  integral  method  over  the 
diffusion  coefficient  calculated  using  the  heat  flux 
equation  (2-58)  or  (2-60).  I  found  that  for  the  surface  120 
cm  of  Site  2  soils,  a  ratio  of  heat  fluxes  HF  ranging  from 
1.0  to  1.8  resulted  in  a  ratio  of  water  vapor  diffusion 
coefficients  DC  ranging  from  1  to  over  30.  This  discrepancy 
is  caused  in  part  by  error  propagation  involved  in  (2-71), 
which  requires  that  two  nearly  equal  numbers  be  subtracted 
(provided  the  temperature  gradient  is  negative)  and  the 
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difference  divided  by 

a  small  num 
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the  ratio  of  the  hydraulic  conductivity  c 
proposed  field  procedure  (Section  3.6. 5. 2 
flux  calculated  from  a  water  balance  equa 
hydraulic  conductivity  calculated  usin 
procedure  but  with  the  water  flux  calcula 
found  that  for  the  upper  120  cm  of  Site 
ranging  from  1  to  10  resulted  in  an  HC  r 
1.1  to  over  1  7. 
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analysis  sim 
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by  Gardner  (1  965)  , 

,  Fliihler  et  al. 
out.  Let  a  measured 
endent  variables  xc 
F  =  f  ( x |  ,  ... 

rror  by  ±Ax(  ,  ±Axz  , 
11  cause  an  error 
the  function  f  in  a 
rder  terms,  as  the 
obtains  (Doebelin, 


1975) 


f  (X<  iAX,  ,  X4±4X4f...,  Xy,±AXy,)  = 

f  (x,  ,  x2,...,  Xy, )  +Ax,  5f/^xf 

+£xzdf/2xz  +  ...  +Axy,^f/5xr)  (3-23) 

The  absolute  error  Ed  is  given  by 
Ed=AF=  jAx,  £f  /£x,  1  +|Axz^f/^xz{  +.  .  .+  lAxyy^f/^Xy,  |  (3-24) 

where  the  absolute  value  signs  are  used  to  avoid  possible 
reduction  of  the  total  error  by  negative  signs  of  partial 
derivatives  associated  with  positive  Ax^'s,  or  vice  versa. 
Thus,  the  actual  error  AF  will  never  exceed  as  long  as 
the  Ax,;  *s  do  not  exceed  their  estimated  values.  The  form  of 
(3-24)  is  useful  because  it  indicates  which  variables  (xj_'s) 
exert  the  strongest  influence  on  the  accuracy  of  the  overall 
result.  The  maximum  error  approach  is  different  from  the 
standard  error  approach  (Beers,  1957;  Lyon,  1970),  in  which 
a  statistical  average  of  the  errors  Ax^  's  liable  to  be 
present  —  usually  the  root-mean-sg ua re  errors  —  is 
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estimated. 

In  order  to  illustrate  the  practical  usefulness  of  this 
technique,  approximate  measures  of  the  uncertainty  in  the 
various  parameters  involved  in  calculating  the  heat  flux 
(3-20)  in  a  general  manner,  as  related  to  the  present  field 
measurements,  are  obtained  using  the  maximum  error  approach 
outlined  above.  The  results  of  this  analysis  are  detailed  in 
Table  3-4.  For  error  analysis  it  is  necessary  to  arrive  at 
an  estimate  of  error  (Ax^ )  for  each  of  the  variables.  These 
error  estimates  could  be  standard  deviations,  95%  confidence 
intervals,  ranges  or  other  measures  of  error.  The  maximum 
absolute  error  estimates  used  here  (column  (4)  of  Table  3-4) 
were  based  on  the  instrumentation  and  experiments  used  for 
this  study.  The  error  estimates  of  specific  heats  cs ,  c0  and 
density  of  organic  matter  ^>o  are  based  on  data  given  by  de 
Vries  (1  963)  .  The  variable  x^  exerting  the  greatest 
influence  on  the  accuracy  of  the  overall  result  is  indicated 
by  an  arrow  in  column  (6)  of  the  table  for  each  of  the  four 
example  calculations  shown  there. 

It  should  be  noted  that  for  this  analysis  no  errors 
were  assigned  to  the  problem  of  measuring  true  field 
conditions  caused  by  instrumentation  emplacement  and  the 
resulting  disturbance  of  the  soil,  field  heterogeneities, 
annual  heat  wave  influences  (Hanks  and  Tanner,  1972),  data 
smoothing,  departure  of  reality  from  theory,  and  possibly 
other  factors.  Because  the  computed  errors  depend  on  the 
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TABLE  3-  4  • 

flaxiaum  error  analysis  for  heat  flux  parameters  at  site  2 


(1)  (2)  (3) 


(<*> 


(5) 


(6) 


(7) 


(8) 


Equation  6  mag-  x^  .lagnitude  of  Approx,  flax.  abs.  (jF/ax^  )  A  x^  Contrib.  to  abs.  Total  abs.  Total  relati?© 

nitude  of  P  error  t  error  of  term  error  error  (in  X) 
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*Data  from  a  depth  of  60  cm  from  Site  2  during  the  August  13-14,1976  aeasurement  period 
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magnitude  of  the  variables  involved  in  the  error  analysis, 
representative  data  from  some  depth  (60  cm  in  the  examples 
shown  in  Table  3-4),  where  the  time  change  in  various 
physical  properties  is  not  very  pronounced,  were  arbitrarily 
chosen.  The  purpose  of  this  exercise  is  not  to  present  a 
precise  and  detailed  error  analysis  but  to  get  a  general 
feeling  of  the  approximate  maximum  errors  involved  in  this 
type  of  field  calculation. 

One  should  not  be  discouraged  by  the  error  analysis 
shown  above  because  any  other  field  method  would  most  likely 
result  in  similar,  if  not  greater,  uncertainties  (Fl’uhler  et 
al,,  1976).  The  maximum  error  analysis  can  be  used  as  a 
criterion  for  selecting  the  best  field  method  from  others 
available.  Such  a  selection  was  done,  for  example,  by  Hanks 
and  Jacobs  (1971)  who  showed  that  the  calorimetric  (  or 
temperature  integral)  method  for  estimating  soil  heat  flux 
resulted  in  a  smaller  relative  error  than  the  heat  flux 
meter  approach  or  the  combination  of  heat  flux  meter  and 


calorimetric  method. 


CHAPTER  4 


MODEL  STUDIES 


4.  1  Introduction 

The  purpose  of  this  chap 
and  experimental  results  of 
the  simultaneous  transport 
subsurface.  This  modeling  e 
reproduce  flow  mechanisms  but 
study  of  these  mechanisms, 
better  understanding  of  the 
establishes  the  framework 
data;  it  is  the  device  to  tes 
theories  through  sensitivit 
means  for  prediction  of  non-i 
thus  for  their  control. 
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4_.  2»_  The  Mathematical  Model 

The  first  step  in  a  quantitative  modeling  effort  of 
simultaneous  heat  and  mass  transfer  in  porous  media  is 
development  of  a  mathematical  model.  Such  a  model  consi 
of  the  system  of  coupled,  non-linear,  second-ord 
parabolic  partial  differential  equations  developed 
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Chapter  2  (2-54  and  2-67)  and  presented  below  in 

one-dimensional  form: 

F  Sv|//£t  =  ^  ( Kvjj  ^/^z)/^z  +^(Dt  3T/az)/^z  -aK/a  z  (4-  la) 

c  avat  =  a(^aT/az)/az  o^^/dz)/^  -c^aT/az  (4-ib) 

where 

F  =  a0/c)^  in  the  unsaturated  zone 
=  Ss  in  the  saturated  zone 

Kvp  =  K  +  K^y 

dt  =  Dr/  +  °Tv 
C  =  0-46xy^w,  +  0.60xo  +  0 

dL  =  l'  k^v  (4- 1c) 

q^/o  =  -K^v|)/dz  -DTdT/dz  +  Ki  (4-1d) 

Because  the  flow  of  heat  and  water  is  treated  as  a 
coupled  phenomenon,  it  is  necessary  to  specify  initial  and 
boundary  conditions  at  the  same  boundaries  for  both 
equations  (4-1a,b).  Therefore,  I  used  the  following  type  of 
initial  and  boundary  conditions  associated  with  (4-1a,b). 

Initial  Conditions  T(z,t)  =  T  (z)  5)  t=0,  z€(0,L)  (4-2a) 

Vp(z,  t)  =  vj^(z)  5)  t=0 ,  Z€(0,L)  (  4-  2b) 


Top  Boundary  Conditions  T(0,t)  =  T(t)  a)  t>0,  z  =  0  (4-3a) 
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qw/^(0,t)  =  ±  E(t)  (4-3b) 

specified  evaporation  or 

infiltration  rate 


Bottom  Boundary  Conditions  ^T/^z|z_L  =0  (4-3c) 

gw/^(Lft)  =  ±  w(t)  (4-3d) 
specified  water  flux  entering 
or  leaving  the  system 

The  solution  of  (4-1a,b) ,  with  the  initial  and  boundary 
conditions  (4-2)  and  (4-3)  respectively,  form  a  mathematical 
basis  for  describing  the  spatial  and  temporal  distribution 
of  temperature  and  pressure  head,  as  well  as  all  relevant 
fluxes  and  transport  coefficients  in  porous  media.  The  form 
of  the  present  mathematical  model  is  more  general  than  the 
one  used  by  Philip  and  de  Vries  (1957)  in  that  the  pressure 
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completely  saturated,  while  the  9-based  approach  is  limited 
to  homogeneous  and  unsaturated  porous  media  (Corey,  1  977)  . 
However,  one  disadvantage  of  the  v|j-based  approach,  which  is 
avoided  in  the  0-based  formulation,  is  the  strong 
non-linearity  of  F  in  (4-1a),  especially  in  the  vicinity  of 
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very  sharp  wetting  fronts. 

4. 3  The  Numerical  Approach 

The  complexity  of  the  system  of  coupled  non-linear 
partial  differential  equations  (2-54)  and  (2-67)  suggests 
that  analytical  solutions  are  not  possible.  Numerical 
techniques  provide  the  most  practical  means  for  obtaining  a 
solution  to  the  (4-1)  system  of  partial  differential 
equations.  Because  the  direction  of  heat  and  moisture  flow 
is  predominantly  unidirectional,  I  chose  to  solve  the 
one-dimensional  form  of  the  coupled  equations.  Besides 
reducing  the  complications  of  the  modeling  effort,  the 
one- dimens ional  form  also  reduces  the  amount  of  experimental 
data  required  and  thus  the  uncertainty  associated  with  their 
calculation.  The  latter  situation  holds  especially  with 
regard  to  determining  transport  coefficients  for 
non-homogeneo us  porous  media.  However,  once  a 
one-dimensional  solution  to  the  relevant  system  of  partial 
differential  equations  is  established,  the  extension  to  two 
or  three  dimensional  situations  should  be  feasible. 


The  system  of  equations  (4-1)  may  be  solved  using  a 
variety  of  numerical  methods,  the  most  common  of  which  are 
the  finite-difference  (Forsythe  and  Wasov,  1960;  Richtmyer 
and  Morton,  1967;  Smith,  1969;  Mitchell,  1969;  von 
Rosenberg,  1969;  Remson  et  al.,  1971)  and  finite  element 
methods  (Zien kiewicz , 1 97 1 ;  Desai  and  Abel,  1972;  Sagerlind, 


1976;  Pinder  and  Gray,  1977).  Although  the  recently  po 
finite  element  methods  are  more  advantageous  for  tw 
three-dimensional  problems,  especially  if  complex  geome 
are  involved,  at  present  they  show  little  or  no  adva 
over  the  better  known  finite  difference  techniques 
non-linear  transient  one- dimens io nal  problems  (Emer 
Carson,  1971;  Mercer  and  Faust,  1976). 


The  finite  difference  approach  was  followed  here, 
derivatives  in  the  coupled  equations  (4-1)  and  i 
corresponding  initial  and  boundary  conditions  (4-2) 
(4-3)  were  replaced  by  finite  difference  approximatio 
discrete  points  of  the  solution  domain.  The  resulting  s 
of  non-linear  algebraic  equations  are  solved  for  ^  an 
The  water  characteristic  input  function  was  utiliz 
obtain  the  corresponding  0  values.  Figure  4-1  sho 
schematic  diagram  of  the  one-dimensional  model  chos 
represent  continuous  water  and  heat  flow  from  the  g 
surface  to  a  point  beneath  the  water  table.  The  subsu 
profile  is  discretized  using  the  equidistant  f 
difference  nodes  shown  in  Fig.  4-1. 


4.2.4  Numerical  Solution 

Although  explicit  methods  for  solving  finite  diffe 
equations  are  simple  and  straightforward,  the  s 
restrictions  on  mesh  size  and  time  increments  that 
sometimes  necessary  make  this  method  rather  unsuitabl 
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Upper  b.c. 


Bottom  b.c. 


T  (t)  =  T*  (t) 
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Figure  4-1.  One-dimensional  coupled  transient  mathematical  model 
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practical  applications.  Thus,  I  have  adopted  an  implicit 
method  which  is  less  restrictive  but  numerically  more 
complicated  because  it  involves  the  solution  of  a  system  of 
algebraic  equations  at  each  time  step.  Specifically,  a  form 
of  the  Cr an k-N icholson  implicit  method  (Richtmyer  and 
Morton,  1967)  was  followed  (Appendix  D) .  This  implicit 
scheme  converges  with  a  smaller  discretization  error, 
0  { (At)  2  +  (Az)  2]  ,  compared  to  the  error  0{(At)  +  (Az)  2]  of  the 
backward-difference  implicit  method.  The  resulting  set  of 
simultaneous  non-linear  algebraic  equations  using  the 
space-time  network  of  Fig.  4-2  is  as  follows 
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Figure  4-2 
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-Al+l-,  +  B -  C^l, 

-  K  T L,  *  Bl  T- 

-  C^L  =  HC 

( 4- 4a) 

-KTli  +  Bc,Ti  -Ct'T^ 

-  W-l  +Bi'+i  * 

Ct'  ^  =  HJ 

( 4- 4b) 

where 

j-'/z 

Ai  ~ 

(4-4c) 

Bi  =  2Ff->r  +  + 

4--1 

(  4-4d) 

Ci  =  K^tti 

(  4-4e) 

r  =  At/(Az)2 

( 4-4F) 

IL  =  D~~Yz 

( 4-4G) 

Bi  =  Di:'^  +  v£>/z. 

(  4-4H) 

*  - 

(4-41) 

Hi.  =  -tfl  + 

^  (4'-'  -4> 

*  04;  -t 

/"> 

+  47^.  <4|  -t/  )  + 

(Ki-.4-4+1i)2^ 

+  (2FlJ’l/Vr) 

(  4-  4  J ) 

Ai1  =  +  cuAi7lAz/ 2 

(4-4K) 

=  2  cfA/r  +^+hl- 

-'A 

-'A. 

(4-4L) 

CJ 

(  4-4M) 

_  j-Vi 

V  =  Duc_/z 

(  4-4N) 

Bl  -  nJ'^  +  DJ"'/l 

Bl  di -L+i/z. 

(  4-40) 

Ci?  =  <l\ 

( 4-4P) 

H*  =  ^J'4  <4.  -Tf)  + 

3H4<*£  -^)  + 

<‘4  -<> 

+  DL^/  (C,1  -^')  -l£*.')/2 

+  (2CJr'/Vr)  rf' 

(4-4Q) 

xi'll  =  <*T  ♦*£  +X-L 

+  XJL  )  /4 

(  4-4R) 

=  <4+xL.  <  )/<* 


(  4-US) 
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J-«  J 
(?L  -YJ/2 


( 4-4T) 


(The  reader  is  referred  to  Appendix  D  for  the  forms  that 
(4-4)  take  at  the  top  and  bottom  boundaries  of  the  solution 
domain.  ) 

The  system  of  equations  (4-4)  may  be  represented  in 
matrix  form  as  follows,  where  brackets  indicate  square 
matrices  and  braces  denote  column  matrices: 

C  G  ]  {v|i}  +  c  F  ]  {T}  =  {H} 

[  2]  m  +  C  D]{T}  =  { H  *  }  (  4-5a) 

or 


Each  of  the  matrices  G,  F,  E,  D  above  represents  a 
tridiagonal,  diagonally  dominant  square  matrix.  Thus  trie 
system  (4-5)  has  the  following  form 
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(  4- 5c) 


This  system  of  algebraic  equations  is  non-linear 
because  the  coefficient  matrices  G,  F,  E,  D  depend  on  values 
of  vj^and  T ^  ,  which  are  unknown  at  the  time  the  solution  to 
(4-4)  or  (4-5)  is  sought.  In  order  to  remove  this 
non-linearity,  an  iteration  process  was  used  in  which 


initial  estimates 
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iW  and 
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Linear 
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At  this  point  let 
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linear 


where 


[A]{x]  =  {b} 

[A]  is  an  nxn  constant  matrix; 


(4-6) 
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McCracken,  1972). 

The  matrix  A  may  be  partitioned  into  three  separate 
matrices 

A  =  L+D+U  (4-7) 

where 

L  is  a  strictly  lower  triangular  matrix  whose  entries 
are  the  entries  of  A  below  the  main  diagonal; 

D  is  a  diagonal  matrix  consisting  of  the  diagonal 
entries  of  A;  and 

U  is  a  strictly  upper  diagonal  matrix  whose  entries 


are  the  entries  of  A  above  the  main  diagonal. 

(4-6)  may  now  be  written  as 

[  L  +  D  +  U  ]{x}  =  {b}  (4-8) 

which  gives  rise  to  the  following  well-known  Jacobi  (4-9), 
Gauss-Seidel  (4-10)  and  successive  over-relaxation  (4-11) 
iterative  schemes  for  the  solution  of  systems  of  linear 
equations  (Varga,  1962) 

xUt  =  D~ 1  [  b-  (L+  U)  x L  ]  (4-9) 

XL+<  =  (L  +  D) [  b-UxL  ]  (4-10) 

XL*'  =  (ajL  +  D)  ~ij>b+  {  (1  -ft))  D-fc>U}  xL  ]  (4-11) 


where  u  is  a  relaxation  factor  which  lies  between  1  and  2, 
and  i  is  an  iteration  index.  Note  that  (4-11)  reduces  to 
(4-10)  when  cJ  =  1.  The  matrices  M  =  -D-1(L+U)  from  (4-9),  M 
=  -(L+D)-1U  from  (4-10)  and  M  =  (ujL  +  D)  ~ 1  {  (1  -w)  D-u>U)  from 

(4-11)  are  defined  as  the  Jacobi,  Gauss-Seidel  and 
successive  over-relaxation  iteration  matrices,  respectively. 
Using  the  above  definitions  of  the  iteration  matrices,  each 
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of  the  previous  three  iterative  methods  (4-9  to  4-11 
the  general  linear  form 

xc**  =  Mxl  +  b* 

A  necessary  and  sufficient  condition  for  convergence 
a  linear  iterative  scheme  is  that  the  spectral  radius 
iteration  matrix  M  —  ^(M)  —  that  is,  the 

eigenvalue  in  absolute  value  of  the  iteration  matrix 
less  than  unity  (Varga,  1962).  By  Ge rschgorin • s 
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(1971)  describe  some  linearization  methods  for  se 
non-linear  problems  in  subsurface  hydrology.  Alt 
predictive-corrective  schemes  (Douglas  and  Jones,  1963 
three-level  schemes  (Lees,  1966;  Mitchell,  1969) 
effective  linearization  methods,  a  Picard-type  iter 
(Remson  et  al. ,  1971;  Whisler  and  Klute,  1965)  was  used 
with  the  objective  of  keeping  the  numerical  solution  s 
as  simple  as  possible.  This  simplicity  criterion  is  als 
reason  for  using  a  uniform  grid  throughout  the  st 
profile.  The  linearization  scheme  used  --  together  wit 
main  steps  of  the  numerical  procedure  --  are  indicate 
the  following  block  diagram  (Fig.  4-3). 
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Figure  4-3 


135 


G 

E 

G 

E 


F 

D 

0 

D 


H 

H  ' 


G 

D 

As  G  an 

Thomas 

1969) 

This  sc 

solutio 

iter ati 

over-re 


d  D 

ar 

algo 

ri 

was 

u 

heme 

P 

n. 

Co 

ons. 

T 

laxa 

ti 

w 

H  -  F  TJ  , 

*  rlr^ 

m  \  i  J 

^  .  J 

11  T  x 

e  diagonally  dominan 

thm  (Carnahan  et  a 

sed  to 

solve  the  s 

roved  to 

be  very 

nver gence 

usually  t 

hus  there 

was  no  nee 

on  technique. 

4. 7  Computer  Simulation  Program 


A  listing  of  the  FORTRAN 
solve  tiie  system  of  equations  (4- 
E.  All  calculations  were  carried 
computer  of  the  University  of  A1 
brief  discussion  of  some  of  the 
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grid  size  Az  =  12  cm.  The  initial  water  table  position  was 
at  300  cm  below  the  surface  node  i  =  2.  The  general 

direction  of  flow  was  upwards.  A  time  step  of  0.5  hr  was 
chosen  by  a  limited  trial  and  error  procedure.  In  most  of 
the  early  runs,  a  time  step  of  .At  =  0.  2857  hr  was  used.  The 
entire  subsurface  profile  was  assumed  to  be  homogeneous  with 
respect  to  the  thermal  conductivity,  which  in  turn  was 
considered  a  function  of  water  content  only.  However,  with 
regard  to  hydraulic  conductivity,  the  saturated  domain  of 
the  subsurface  profile  was  treated  as  a  heterogeneous  medium 


Kc 


k'l-M 
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L-H/z 
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homogeneous.  I 

a  flow 

theory  (Scott 

ivit  y 

between  adjacen 
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.  4-4). 

2  A  z 

{(Az/Kt+|  )  +  (Az/KL)} 


(4-14) 


Figure  4-4 


The  saturated  zone  was  divided  into  three  layers  (Fig.  4-1) 
with  field- measured  hydraulic  conductivities  of  K  =  1.256 
cm/hr,  for  the  first  layer;  K  =  0.687  cm/hr  for  the  second, 
and  K  =0.054  cm/hr  for  the  bottom  layer,  which  indicated  the 
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presence  of  a  loamy  lens  at  that  position.  Hydraulic 
conductivities  were  corrected  for  the  effects  of  temperature 
on  density  and  viscosity  using  the  following  relationship 

K<U  =  ‘WV  (4-15) 

where  the  reference  temperature  T  was  IOC  for  the 

saturated  zone  and  =  20C  for  the  unsaturated  zone. 

The  porous  medium  was  assumed  to  possess  a  constant 
porosity  value  =  0.41  and  elastic  storage  properties  that 
can  be  effectively  represented  by  the  specific  storage  Ss  . 
The  specific  storage  was  assigned  a  constant  value  (Ss  = 
2x10~3  cm-1)  for  the  saturated  region  —  which  is  an 
acceptable  approximation  when  the  changes  in  water  table 
elevation  occur  slowly  (McWhorter  and  Sunada,  1977)  —  and 

zero  for  the  unsaturated  zone.  In  effect,  the  unsaturated 
porous  medium  was  assumed  to  be  incompressible. 


The  top  boundary  conditions  consist  of  a  soil  surface 
temperature  specification  for  the  heat  flux  equation  and  a 
prescribed  water  flux  (infiltration  or  evaporation)  for  the 
pressure  head  equation.  The  surface  temperature  may  be 
specified  as  a  sinusoidal  wave  or  as  piecewise  functions 


joining 

meas  ured 

data.  In  the  case  of 

evaporation  condition 

the  user 

does  not 

have  to  specify 

an 

evaporation  flux 

relation 

because 

it  is  calculated 

in 

the  program.  The 

evaporation  flux  is  calculated  by  using  the  water  balance 
equation  (3-21)  with  the  input  flux  calculated  using  the 
non-isother ma 1  flux  equation  (2-38),  and  utilizing  the  water 
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contents  of  the  current  and  previous  time  steps  to  calculate 
the  time  change  of  moisture  content  for  the  top  mesh  layer. 

The  bottom  boundary  water  flux  was  given  a  range  of 
expected  values  based  on  direct  or  indirect  observations  of 
the  hydraulic  gradient  and  conductivity  of  the  bottom 
material.  The  value  assigned  to  the  bottom  boundary  flux  was 
w  =  0.07  cm/hr  downward.  The  sign  convention  followed  was 

negative  for  upward  flow  and  positive  for  downward  flow.  The 
bottom  Doundary  for  the  heat  flux  equation  was  a  zero 
temperature  gradient  specif ication.  This  gradient  was 
approximated  by  a  forward  difference  scheme,  in  contrast  to 
the  central  difference  schemes  used  generally  throughout  the 
computer  program  to  approximate  spatial  derivatives. 

As  a  convergence  criterion,  a  tolerance  of  6=0.  25  or 

0.3C  was  used  for  the  computer  runs  with  a  At  =  0.2857  hr 

and  €  =  0.  5C  for  those  with  a  At  =  0.5  hr  between 

temperatures  calculated  during  two  successive  time  steps, 

excluding  the  surface  temperatures.  Therefore,  if  the 

)  JH 

condition  |T^  T^  |  <6  was  satisfied,  the  computed 

solution  was  accepted  as  valid;  otherwise,  an  iteration 
cycle  was  activated  (Fig.  4-3).  A  similar  convergence 
criterion  could  have  been  established  for  the  vp-values. 
However,  judging  from  the  results  obtained  using  only  the 
temperature  convergence  criterion,  such  a  condition  was  not 
considered  necessary  because  the  pressure  head  values  may 
change  considerably  in  a  time  step,  especially  during 
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infiltration 

into 

a  dry  soil. 

The 

numerical 

sche  me 

I  used 

proved  to 

be 

very  fast 

in 

con  ve  rging 

to  the 

above 

require  ment  s 

with 

only  one  or 

t  wo 

iterations. 

If  the 

n  umber 

of  iterations  exceeded  a  prescribed  number  KW(g(  without 
convergence,  the  procedure  was  stopped  with  appropriate 
messages  printed  (Fig-  4-3).  The  computer  simulation  was 
also  stopped  and  a  warning  message  printed  if  the  value  of 
the  bottom  node  became  negative. 


Because  of  the  banded  nature  of  the  square  matrices  D, 
E,  F,  G  (4-5c)  with  one  upper  and  one  lower  diagonal,  these 
matrices  were  stored  row-wise  (oand  storage  mode)  so  that 


the 

zero 

elements  were 

compressed 

out 

of  the  matrices.  Thus, 

the 

above 

nxn  matrices 

were  reduc 

ed  t 

o  nx 3  matrices 

with  the 

diagona 1 

elements  falling  into 

the 

second  column. 

Equation 

(4-5)  shows  an  example  of  the  band  storage  mode  used  to 
minimize  computer  memory  requirements. 
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The  simulation  program  can  be  controlled  by  the  user  in 
following  ways.  By  setting  the  computer  input  parameter 


the 


' 
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KODE1  equal  to  1  or  0,  the  bottom  boundary  condition  can  be 
specified  as  recharge  or  discharge.  By  setting  KOBE 5  equal 
to  0  or  1 ,  the  top  boundary  can  be  specified  as  evaporation 
or  infiltration.  In  a  similar  way,  by  putting  KODE4  equal  to 
0  or  1,  the  user  may  specify  a  problem  either  in  an 
uncoupled  (isothermal)  or  coupled  (non-isot hermal)  transient 
(unsteady)  form.  KODE2  =  i,  where  i  =  1,  2,  ...,  NJ  ( NJ  = 
number  of  time  steps)  specifies  the  number  of  time  steps 
that  may  be  skipped  before  the  result  is  printed. 

At  the  beginning  of  the  programming  effort,  balance 
equations  were  used  to  derive  the  hydraulic  conductivity 
functions.  The  water  flux  was  calculated  from  the  measured 
time  change  of  water  content  and  the  known  water  flux  at  the 
bottom  of  the  top  mesh  layer.  A  similar  balance  equation  was 
also  used  to  calculate  the  heat  flux  from  the  measured  time 
change  of  the  soil  heat  content  with  the  bottom  heat  flux 
calculated  from  the  thermal  conductivity  and  temperature 
gradient.  However,  these  schemes,  especially  the  latter  one, 
resulted  in  numerical  oscillations  and  inaccuracies. 
Therefore,  the  field  approach  proposed  in  Chapter  2  was  not 
implemented  in  the  present  computer  program. 

As  a  consequence  of  the  above  development,  the 
transport  coefficients  involved  in  the  system  of  equations 
(4-1)  were  estimated  analytically  as  shown  in  Chapter  2.  In 
order  to  increase  the  accuracy  of  such  calculations,  I  used 
double  precision  arithmetic  throughout  the  program.  Because 


141 

of  uncertainties  surrounding  the  D-p^  coefficient  for 
saturated  conditions  mentioned  in  Chapter  2  (p.  29)  ,  and  in 
order  to  simplify  the  numerical  calculations  somewhat,  the 
thermal  diffusion  coefficient  D-p  was  assumed  equal  to  zero 
under  saturated  conditions.  This  assumption  does  not  seem 
unreasonable  in  view  of  the  small  temperature  gradients 
observed  in  the  saturated  zone  under  study  and  of  its  small 
thickness;  also  the  limited  number  of  relevant  data  in  the 
saturated  zone  severely  limited  direct  use  of  the  quantity 
■^/  c)  T  (2-31)  by  introducing  considerable  errors. 

The  ratio  £  of  the  average  temperature  gradient  in 
air-filled  pores  to  the  overall  temperature  gradient,  which 
appears  in  the  expression  of  D-pv  coefficient  (2-18),  has  to 
be  specified.  The  following  values  were  adopted  from  data 
provided  by  Philip  and  de  Vries  (1957)  and  Rose  (1968a,  b) . 
For  T  <  25C,  <£  =  1.8,  and  for  T  >  25C,  1.3.  It  should  be 

noted  that  a  conservative  point  of  view  was  adopted  in 
choosing  the  ^  values  and  in  incorporating  Rose’s  £ 
coefficient  —  the  liquid-assisted  vapor  enhancement  factor 

—  into  the  pore  geometry  factor  f.  The  computer  program 
provides  for  both  vj^-based  and  6-based  equations  discussed  in 
previous  chapters. 

Experimentally  determined  input  data  to  the  simulation 
model  include: 

—  initial  distribution  of  temperature  and  pressure  head 
with  depth  (usually  supplied 


as  graphically  smoothed 
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curves)  ; 

—  time  distribution  of  surface  t 
in  Chapter  3,  the  measured  s 
of  less  than  1  cm  were  taken  a 

—  characteristic  function  (s)  -- 

—  hydraulic  conductivity  funct 
function  may  be  optional  bee 
use  the  water  balance  equation 
into  the  model) ; 

—  thermal  c onductivity- water  c 
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physical  parameters 

for  the 

system,  sue 

h 

as  the 

section 

length  studied,  the  gas 

constant,  th 

e 

specific 

storage 

value,  and  others; 

numerical  informat 

io 

n,  such 

as  the 

maximum  simulation 

time , 

the  size  and 

n 

umber  of 

time  and 

space  increments. 

the 

maximum  numb 

er 

of  allowable 

iterations,  and 

others; 

a  list 

of 

top 

boundary 

time-temperature  distribution  and  initial  depth  distribution 
of  6  and  T;  the  calculated  temperature,  pressure  head  and 
moisture  distribution,  as  well  as  all  transport  coefficients 
and  heat  and  water  fluxes. 

As  the  above  discussion  may  suggest,  this  simulation 
model  is  general  and  flexible.  It  can  handle  any  type  of 
boundary  conditions,  heterogeneities,  saturated- unsaturated 
regions,  isothermal  and  non-isothermal  environments, 
transient  and  steady-state  systems. 

4._8  Program  Verification 

The  overall  error  of  a  simulation  result  may  be 
considered  as  the  sum  of  errors  introduced  by  the  numerical 
finite  difference  analog  and/or  errors  introduced  by 
possible  inadequacies  of  the  theoretical  framework  in 
describing  the  particular  phenomena.  Although  the  validity 
of  the  uncoupled  (isothermal)  equations  is  supported  by  a 
considerable  body  of  literature,  the  correspondence  of  the 
coupled  non-isothermal  equations  to  the  physical  processes 


has  not  been  well  established.  Such  uncertainties  are  mainly 


due  to  the  relative  scarcity  of  stud 
validity  of  the  non- isothermal  coupled 
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several  aspects  of  such  systems.  However 
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algorithm  used  in  the  present  study 
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A  limited  number  of  analytical  solutions  are  avai 
for  the  complete  heat  transport  equation.  However,  no 
analytical  solutions  for  the  moisture  equation  are  avai 
(Guymon  and  Luthin,  1974).  Kirkham  and  Powers  ( 
summarize  some  of  Philip's  (1969)  approximate  sol 
techniques  of  uncoupled  moisture  problems.  Therefor 
order  to  check  the  validity  of  the  numerical  method. 
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transient  uncoupled  heat  equation  with  constant  coefficients 
was  solved  using  the  same  computer  model.  The  results  were 
compared  with  the  analytical  solutions  described,  for 
example,  in  Carslaw  and  Jaeger  (1959)  or  in  Ingersoll  et 
al.,  (1954).  Specifically,  the  following  heat  conduction 
problem  was  considered: 

C  3T/3t  =  3(3  ST/3z)/dz  =3(^2T/3z2)  (4  -1 7a) 

dT/3t  =  (3/C)  32T/3z2  =  OC  32T/^z2  (4-1 7b) 

where  o(  is  the  thermal  diffusivity  (cm2/hr)  with  the 
following  boundary  and  initial  conditions: 

T  =  OC  a)  z  =  0 


T  =  OC  a)  z  =  L 


T  =  T0  =  IOC  a)  t  =  0 


(4-1  7c) 


The  analytical  solution  to  the  above  problem  is 
(Carslaw  and  Jaeger,  1  959): 

(U-18a) 

>1=0  1— 
or 

-q/iVt/r1- 


3jiz.  , 

Sir)  4.  J_  Q. 

L  5 


sin  f'T — 


.}  (4-1 


-18b) 


T=(4T;A){e  '  siv,^  +  ±e  L.  ■  5  -  -  -  L. 

Using  the  following  representative  input  data,  °<  =  J3 /C  = 
8/0.5  =  16  cm2/hr  and  L  =  564  cm,  one  finds  that  for  t  =  24 
hr  and  z  =36  cm,  for  example,  the  analytical  solution  -- 
using  the  terras  shown  in  (4  —18b)  --  is  T  (z*3e,t  =  Z4')  8.07 2, 

while  the  numerical  solution  obtained  using  the  present 
program  with  a  time  step  At  =  0.2857  hr  was  T  (z  =  36_,  = 

8.05C.  The  agreement  between  the  two  is,  therefore, 
excellent  (relative  error  ->0.25%)  suggesting  that  the 
uncoupled  transport  equations  are  solved  correctly.  Figure 


. 
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4-5  shows  the  initial  and  simulated  transient 
temperature-depth  distribution  for  the  conduction  problem 
(4-17)  considered  here. 

The  simulation  results  were  also  compared  with  observed 
field  data.  Figures  4-6  and  4-7  compare  field  and  simulated 
results.  In  all  those  figures  the  input  depth  distribution 
of  temperature  or  moisture  in  the  early  afternoon  of  August 
13th,  1976  at  Site  2  are  shown  with  the  continuous  curves. 
If  one  keeps  in  mind  that  in  the  field  the  unsaturated 
domain  included  some  discontinuous  clay  loam  lenses  but  that 
the  simulation  model  assumed  a  homogeneous  unsaturated  zone, 
the  simulated  results  are  in  satisfactory  agreement  with  the 
observed  data. 

The  present  modeling  study  does  not  represent  an 
exhaustive  numerical  treatment.  Because  of  the  inherent 
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the  accuracy  of  the  finite  difference  approximation  is  a 
function  of  the  fineness  of  the  time  and  spatial  grids  used 
in  the  solution,  subject  to  computational  limitations,  such 
as  truncation  and  round-off  errors.  However,  in  order  to  use 

rather  coarse 


the  available  computer  funds  most  efficiently. 
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4-5.  Time  distribution  of  temperature  for  the 
conduction  problem 
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.  Comparison  of  observed  and  simulated 
temperature  distributions 
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Figure  4-7.  Comparison  of  observed  and  simulated  moisture  distributions 
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grid  systems  were  used  to  reduce  the  storage  and  time 
requirements.  A  computer  simulation  run  of  168  hr  (one  week:) 
with  a  time  step  of  0.5  hr  used  approximately  25  sec  of  CPU 
time  in  the  University’s  fast  AMDAHL  computer;  while,  the  24 
hr  runs  with  a  time  step  of  0.2875  hr  (with  some  minor  later 
modifications)  used  approximately  12  to  13  sec. 
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divergence  of  the  solution.  For  better  and  faster  results, 
the  input  data  sets  should  be  representative,  continuously 
varying  and  smoothed. 

4.9  Results  and  Discussion 

A  limited  number  of  simulation  problems  were  solved 
using  the  computer  simulation  program  described  above.  In 
every  case,  both  the  isothermal  (uncoupled)  and 
non-isother ma 1  (coupled)  models  were  used  in  order  to 
evaluate  possible  coupling  effects.  Three  different  initial 
conditions  were  used  in  studying  the  drying  or  evaporation 
problem.  These  initial  conditions  are  shown  in  Fig.  4-8;  the 
corresponding  computer  runs  are  referred  to  as  Run  1,  Run  2 
and  Run  3.  The  upper  boundary  temperature-time  distributions 
used  for  these  simulations  are  shown  in  Fig.  4-9.  Table  4-1 
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summarizes  the  conditions  for  each  run.  A  list  of  some 
physical  parameters  and  numerical  information  for  these 
simulation  runs  are  shown  in  Table  4-2. 

Using  the  initial  conditions  for  Run  1  (wet  condition 
with  field  observed  but  smoothed  data)  and  for  Run  2  (dry 
condition)  indicated  above,  no  appreciable  changes  in  the 
temperature  field  occurred  by  coupling  or  uncoupling  the 
system  of  equations  (4-1)  for  short  time  simulations. 
However,  significant  differences  were  observed  in  the 
evaporation  and  moisture  fluxes  as  indicated  in  Tables  4-3 
and  4-4. 

It  should  be  noted  in  those  tables  that  the  evaporation 
and  moisture  fluxes  computed  from  the  uncoupled  equations 
are  higher  than  the  ones  computed  from  the  coupled 
equations,  and  that  the  difference  is  greater  when  the  water 
contents  are  at  intermediate  levels  (Run  1).  This  situation 
holds  true  as  long  as  the  normal  temperature  gradients 
prevail  in  the  soil;  that  is,  when  the  higher  temperatures 
are  at  the  top  and  decrease  with  depth,  in  opposition  to  the 
potential  gradients  observed  at  the  field  sites.  Thus,  a 
small  flow  is  set  up  opposite  the  potential  gradient  to  slow 
the  drying  rate.  However,  during  the  diurnal  cycle,  when 
temperature  gradients  are  reversed,  the  opposite  holds  true, 
as  shown  in  Tables  4-5  and  4-6. 

Figures  4-10  and  4-11  indicate  a  much  faster  decrease 
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Figure  4-8.  Initial  conditions  for  case  simulation  runs 
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Figure  4-9.  Upper  boundary  time  distribution  of  temperature 
for  case  simulation  runs 
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equations  for  Runs  1  and  2.  Soil  temperatures  decreasing  with  depth. 
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(p.  160  and  thereafter).  These  figures,  in  conjunction  with 
the  moisture  flux  equation  (2-38) l,  show  that  the  drier  the 
soil  becomes  the  more  dominant  the  second  term  B  in  (4-19) 
becomes.  Therefore,  under  very  dry  conditions  and  with  the 
temperature  decreasing  with  depth,  the  positive  (downward) 
term  B  in  (4-19)  exceeds  the  negative  (upward)  term  A  in 
magnitude.  As  a  consequence,  evaporation  becomes  positive; 
that  is,  a  condensation  effect  results.  This  situation  was 
simulated  by  imposing  the  initial  and  top  boundary 
conditions  of  Run  3  (Figs.  4-8  and  4-9  )  consisting  of  a 
highly  negative  ^  distribution  with  depth,  coupled  with  high 
surface  temperatures  and  a  constant  initial  temperature 
gradient  of  0.55  C/cm  throughout.  A  sample  of  simulated 
results  are  shown  in  Tables  4-7  and  4-8. 


The  transport  coefficients  used  in 

th 

is 

model 

dif 

from  the  familiar 

form  of 

Philip 

an  d 

de 

Vries  * 

(19 

coefficients  in  that 

they  are 

derived  f 

ro  m 

the 

more 

ge  ne 

'j'-based  equations. 

As  can 

be  seen 

i  n 

Fig 

.  4- 1 

0, 

model-calculated  K^v  values  (2-17)  are  extremely  small. 
However,  in  the  range  of  water  contents  where  the  soil 
relative  humidity  h  is  below  100%,  the  considerable  increase 
of  the  pressure  head  gradient  will  cause  the  product 

(Kvv,  ly/b  z )  appearing  in  (4-1)  to  rise  to  a  maximum  value 


iln  one-dimensional  form,  this  equation  becomes 
qw/e  =  '  K0.  3j'/^zj  -  Dt  3T/jz  +  Ki 

^  V  ^  1  V  J  ^  "V  * 

A  B  C 


(4-19) 
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Figure  4-10.  Total  water  (hydraulic)  conductivity  function  (e) 
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Figure  4-11.  Thermal  diffusivity  function  (0) 
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before  eventually  falling  to  zero  as  j)  ,  and  therefore  h, 

approach  zero.  Notice  in  Fig.  4-10  that  dominates  over  K 

at  very  low  moisture  contents. 

The  thermal  diffusion  coefficients  calculated  from  the 
model  are  shown  in  Fig.  4-11.  The  D-j-y  coefficient  increases 
continuously  with  decreasing  moisture  content  and  becomes 
the  dominant  thermal  diffusion  coefficient  at  low  moisture 
contents,  as  can  be  seen  in  the  figure.  The  break  in  slope 

in  the  D-^  -  &  curve  (marked  by  an  arrow  in  the  figure)  is  due 

to  the  break-up  of  the  ^-values  used  in  this  program  (p. 

141).  As  mentioned  before,  there  is  some  uncertainty  in  the 
determination  of  D ^  coefficient,  as  proposed  by  the  Philip 
and  de  Vries  surface  tension  model  (p.  29).  Values  of  the 

coefficient  computed  from  (2-37)  were  compared  with  the 
values  obtained  by  the  Philip  and  de  Vries  model  (2-35)  and 
were  found  37%  to  41%  larger  than  the  latter,  as  shown  in 
Table  4-9.  My  approximation  seems  to  be  a  better  one  in  view 
of  the  fact  that  experimentally  measured  values  of  D  were 
up  to  twice  (Wilkinson  and  Klute,  1962)  or  more  as  large 
(Jury  and  Miller,  1974)  as  those  computed  by  the  Philip  and 
de  Vries  model.  The  D-^  coefficient  (Fig.  4-11)  increases 
with  increasing  water  content  before  starting  to  rapidly 
decrease  near  saturation. 


An  explanation  of  this  coefficient  pattern  is  as 
follows  (Joshua  and  de  Jong,  1973).  Most  of  the  observed 
deviations  of  water  flux  in  porous  media  from  Darcy's  law 
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are  attributed  to  interactions  between  the  water  and  the 
porous  medium  (S war tzendruber ,  1969).  This  interaction  is 
proportional  to  the  area  of  contact  between  the  porous 
medium  and  liquid  water  and  increases  as  the  medium  pore 
sizes  decrease.  Hence,  as  the  water  content  in  the  soil 
decreases,  the  proportion  of  water  flow  that  occurs  in  the 
smaller  pores  increases  and  thus  soil-water  interaction 
increases.  In  general,  coupling  phenomena  are  also 
pronounced  when  interaction  between  the  fluid  and  the  porous 
medium  occurs.  This  phenomenon  can  be  demonstrated  by  the 
increase  in  thermally  induced  water  flow  which  occurs  with 
increasing  soil  dryness,  a  case  presented  in  this  thesis. 
When  soil-water  flow  takes  place  at  high  water  contents,  a 
large  proportion  of  this  flow  occurs  through  the  larger 
pores,  thereby  minimizing  the  effect  of  soil-water 
interaction.  Therefore,  the  decrease  in  Dn  at  high  water 
contents  may  be  attributed  to  low  soil-water  interaction. 
However,  at  very  low  water  contents,  when  vapor  flow  becomes 
of  importance,  the  largest  proportion  of  vapor  flow  would 
occur  in  the  large  air-filled  pores;  thus,  little  coupling 
between  heat  and  moisture  can  be  expected. 


Dirksen  (1969)  tried  to  apply  the  above  observations  to 
saturated  clay  suspensions.  He  reasoned  that  by  decreasing 
the  water  content  of  tnese  saturated  clay  suspensions  by 
increasing  the  compaction  pressure,  the  temperature 
gradients  would  become  more  effective  than  hydraulic 
gradients  in  causing  water  to  move  within  compacted 
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sediments.  Although  he  observed  an  increase  in  thermally 
induced  flow  with  increasing  compaction  pressure,  the  amount 
of  such  flow  in  his  experiments  was  generally  less  than  that 
reported  for  liquid  water  in  unsaturated  soils  at  comparable 
water  contents.  This  result  was  attributed  to  the  influence 
of  the  vapor  phase  which  is  absent  undsr  saturated 
conditions. 

The  calculated  secondary  maximum  in  the  Dr^  - Q 
relationship  (Fig.  4-11)  at  the  low  water  content  range  is 
due  to  the  faster  relative  increase  of  | In  hj  with 
decreasing  water  content  compared  to  a  slower  relative 
decrease  of  K  in  that  range  of  water  contents  (2-37).  Data 
to  check  the  physical  significance  of  this  observation  are 
not  available. 

The  transport  coefficient  Dj_  (4-1c)  appearing  in  (4-1b) 
is  usually  extremely  small  in  magnitude.  As  can  be  seen  in 
Fig.  4-12,  it  increases  significantly  at  low  moisture 
contents.  This  increase  --  combined  with  a  large  increase  in 


the  pressure  head  gradient 


will  cause  this  latent 


heat  flux  component  (D^  H'/a*)  appearing  in  (4-1b)  to  reach 
a  maximum  value  before  eventually  falling  to  zero  as  $ 
approaches  zero. 

The  temperature  field  in  the  simulation  runs  presented 
here  was  not  appreciably  influenced  by  coupling  or 
uncoupling  because  in  the  range  of  water  contents  over  which 
most  of  the  simulations  were  run,  that  is,  above  5%  by 
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volume,  the  conductive  term  of  the  heat  flux  equation  ( 
was  dominant  over  the  latent  and  convective  terms  o 
total  heat  flux.  Under  very  dry  conditions,  however, 
latent  heat  term  of  the  heat  flux  equation  (2-58)  be 
increasingly  important.  At  the  same  time,  the  condu 
term  decreases  rapidly.  It  should  be  noted  that 
convective  term  of  the  heat  flux  equation  (the  third  te 
2-58)  --  which  is  usually  neglected  in  non-isoth 

studies  —  becomes  increasingly  important  at  higher 
contents,  as  shown  in  Table  4-10. 
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As  can  be  seen  in  Figs.  4-10,  4-11  and  4-12 

above,  the  computed  transport  coefficients  are  extr 
variable  with  changing  water  content.  This  non-linearit 
the  system  (4-1)  prevents  characterization  by  con 
coefficients.  It  should  also  be  noticed  that  when 
volumetric  water  content  becomes  less  than  about 
thermally  driven  moisture  is  larger  than  the  pote 
gradient  flow  (Fig. 4-10)  for  the  Taber  sites. 

The  hydraulic  conducti vity-wat er  content  relatio 
used  in  any  water  flow  model  exerts  a  significant  infl 
on  the  simulation  results  because  the  calculated 
fluxes  and  several  transport  coefficients  depend  on 
relationship.  For  the  present  study,  conductivity  func 
for  two  different  depths  of  Site  2  were  experimen 
determined  (Figs. 3-25,  C-14).  Both  hydraulic  conduct 

functions  were  input  to  the  computer  simulation  mode 
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Figure  4-12.  (e)  coefficient  function 
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determine  which  one  resulted  in  a  better  agreement  between 
observed  and  simulated  results.  I  found  that  the  hydraulic 
conductivity  function  of  the  surface  -«20  cm  of  soil  resulted 
in  a  comparatively  better  match  between  observed  and 
simulated  data  and  was  thus  considered  more  representative 
of  the  whole  unsaturated  profile. 
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also  observed  with  the  heat  flux  distribution  calculated  by 
the  model,  as  shown  in  Fig.  4-14.  The  time  distribution  of 
the  pressure  head  in  the  vicinity  of  the  water  table  for  the 
same  drying  conditions  (Run  2)  is  shown  in  Fig.  4-15,  where 
the  gradual  decrease  in  the  water  table  depth  with  time  is 
ap  parent. 

The  evaporation  rate-time  distribution  under  a 
diurnally  varying  environment  resulted  in  a  series  of 
diurnal  evaporation  phases,  as  shown  in  Fig.  4-16,  which 
tended  to  damp  down  in  time  as  the  surface  zone  desiccation 
progressed  to  greater  depths.  Such  a  pattern  was  actually 
measured  in  the  field  by  Kimball  and  Jackson  (1971).  In 
simulating  this  evaporation  pattern,  a  sinusoidal  surface 
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Figure  4-13.  Time  distribution  of  evaporation  flux  and 
temperatures  at  various  depths 


cal/cm2/hr 
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Figure  4-14.  Time  distribution  of  heat  fluxes  at  various  depths 


►  cm 
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Figure  4-15.  Time  distribution  of  pressure  head  in  the 
vicinity  of  the  water  table 
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temperature  distribution  similar  to  the  one  shown  in 
4-13  and  the  initial  conditions  for  Run  2  (Fig.  4-8) 
applied.  It  should  be  noticed  that  after  the  third  d 
evaporative  drying,  a  stable  evaporation  rate-time 
pattern  was  observed.  The  time  distribution  of  the 
flux  at  12  cm  depth  is  also  indicated  in  Figure  4-16. 
distribution  exhibits  a  similar  diurnal  cycle  but 
smaller  fluctuations  decreasing  progressively  in  the  d 
layers. 

Although  special  emphasis  was  placed  on  the  d 
problem,  the  infiltration  problem  could  also  be  handle 
the  present  simulation  model.  An  example  of  a  we 
situation  is  given  in  Fig.  4-17,  where  the  time- 
pressure  head  profiles  under  a  constant  infiltration  ra 
0.09  cm/hr  are  shown. 
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Figure  4-16.  Cyclic  evaporation  distribution 
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cm  H20  -* - 0 


Figure  4-17.  Pressure  head  distribution  during  infiltration 


CHAPTER  5 


CONCLUDING  COMMENTS  AND  RECOMMENDATIONS 


The  present  study  has  demo 
applying  and  modeling  the  non-iso 
field  situation.  Comprehensive 
collected  here  are  generally  unav 
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This  work  showed  that  the  mechanistic  Philip  and  de 
Vries-type  of  non-isothermal  theory  can  be  modified  and 
extended  to  saturated  conditions.  These  developments 
involved  the  introduction  of  a  new  approach  for  the 
calculation  of  liquid  thermal  diffusion  coefficients  based 
on  local  vapor-liquid  equilibration  in  porous  media  (2-37) ; 
also  of  procedures  for  calculating  the  various  transport 
coefficients  based  on  calculating  ratios  of  these 
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coefficients  and  combining  these  ratios  with  a  detailed  heat 
flux  equation  (Section  2.8).  Several  problems  involving  the 
application  of  the  proposed  field  methodology  became 
apparent,  such  as  the  inadequacies  of  the  thermocouple 
psychro mete rs  used  in  the  field,  the  insensitivity  of  some 
instruments  to  small  changes  in  field  parameters,  the  very 
small  magnitude  of  various  transport  coefficients  involved 
in  non-isot he rmal  flow  studies,  and  others  (Section  3.8). 


A  comparative  study  of  the  mechanistic  and  irrever 
thermodynamics  approach  to  coupled  problems  indicated 
the  latter  has  not  yet  advanced  beyond  carefully  contr 
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accurate  measurements  of  heat  and  water  fluxes  are  esse 
to  test  the  validity  of  non-isother ma 1  equations. 

This  work  showed  how  field  equipment  can  be  improv 
modified  to  meet  the  particular  needs  of  the  researcher 
example,  in  order  to  meet  the  challenge  of  the 
experimentation,  novel  procedures  in  design  and  constru 
of  field  equipment  were  developed  with  respect  to  th 
conductivity  probes,  piezometers  and  the  va 

combinations  of  data  recording  devices.  Thus,  a  new  d 
for  installing  and  recovering  the  thermal  conduct 
probes  was  implemented;  also,  minor  modifications  in 
construction  of  thermal  conductivity  probes  were  ca 
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This  study  was  the  first  to  demon stra 

possibilities  for  using  a  unified  unsa tura ted-sa 
approach  to  non-isotherma 1  flow  conditions.  As  a  resu 
scope  or  scale  of  applicability  of  the 

non-isother mal  transport  parameters  and  equations 
were  previously  restricted  to  the  uppermost  layer 
ground  surface  or  small  lab  columns  —  was  exten 
accommodate  such  an  approach. 

My  research  showed  that  it  is  possible  to  mo 
simultaneous  flow  of  water  and  heat  in  unsa tura ted-sa 
domains  based  on  an  adaptation  of  the  Philip  and  de 
model.  I  developed  and  tested  this  modification 
measured  field  and  laboratory  data.  I  found  the  ag 
between  the  transient  numerical  solution  of  the 
equations  under  varying  top  boundary  conditions  a 
field  observed  data  satisfactory.  The  present  model  p 
soil  temperatures,  pressure  heads  and  water  contents 
as  all  relevant  transport  coefficients  and  fluxes 
general  enough  to  be  applied  to  related  proble 
simultaneous  heat  and  water  transfer  where  a  pre 
capability  is  needed. 

The  influence  of  coupling  on  diurnal  subsurface 
and  heat  transfer  can  be  determined  by  compari 
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non-isother mal  and  the  isothermal  solutions  of  the  wate 
heat  flow  equations  using  the  appropriate  options  built 
the  present  simulation  model.  The  results  confirmed 
coupling  did  not  noticeably  influence  the  temperature 
but  did  affect  the  evaporation  and  subsurface  moi 
fluxes.  Thus,  depending  on  the  direction  of  the  temper 
and  pressure  head  gradients,  the  subsurface  moisture  f 
might  have  been  underestimated  or  overesti 
significantly  (up  to  -*40%  for  cases  studied  here)  usin 
isothermal  model.  The  analysis  of  transport  coefficient 
the  study  area  indicated  that  thermally  induced  moi 
movement  becomes  dominant  at  very  low  moisture  contents 

The  non-isothermal  theory  can  be  used  to  calc 
evaporation  from  the  soil  surface  based  on  moistur 
temperature  profiles.  Comparison  of  evaporation 
calculated  using  the  coupled  and  uncoupled  equa 
indicated  that  non-isothermal  evaporation  rates  were  h 
or  lower  than  the  isothermal  ones  depending  on  the  dire 
of  the  temperature  and  pressure  head  gradients, 
magnitude  of  such  differences  was  up  to  a  factor  of  two 
the  cases  studied  here,  depending  on  the  soil  moi 
levels.  Thus,  for  example,  hourly  evaporation  esti 
using  an  isothermal  model  could  be  overestimated  by  a  f 
of  two. 

A  cyclic  evaporation  pattern  was  evident  under  di 
conditions.  As  a  consequence,  the  three  stages  of 
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evaporation  process  (Lemon,  1956)  under  constant 
evaporativity  (  or  evaporative  demand  of  the  atmosphere) 
lose  their  distinctness.  At  the  end  of  a  week-long 
simulation,  such  cyclic  evaporation-time  waves  were  still 
apparent. 

Because  of  the  predominance  of  the  conductive  heat 
transfer  over  all  other  modes  of  heat  transport  in  the  cases 
studied  here  and  the  generally  small  magnitude  of  the 
thermal  diffusion  coefficient  compared  to  the  hydraulic 
conductivity  coefficient  under  wet  soil  conditions,  the 
isothermal  model  may  be  used  adequately  to  predict  pressure 
head  and  temperature  profiles  on  a  short-term  basis  (greater 
than  24  hr).  However,  for  long  term  simulations  under 
extended  drying  conditions,  the  coupled  model  should  be 
preferred  in  studying  the  long  range  seasonal  influences  on 
water  transport,  especially  in  view  of  the  fact  that 
measured  values  of  the  D_^  coefficient  are  much  larger  than 
model  predicted  ones  (Jury  and  Miller,  1974). 


Much  remains  to  be  done  in  order  to  consolidate  and 
expand  the  applicablity  of  coupled  non-isothsrmal  models. 
Further  generalization  of  results  obtained  here  will  require 
the  study  of  environments  differing  with  respect  to  the 
porous  materials  and  the  ranges  of  pressure  heads.  More 
experimental  information  is  needed  on  the  thermal  dependence 
of  the  pressure  potential,  particularly  under  saturated 
conditions.  Quantification  of  such  information  would  assist 
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Figure  A-l.  Time-depth  distribution  of  temperature  at  Site  1 
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Figure  A-2.  Potential  evaporation  fluxes,  air  temperatures  and 

wind  speeds  at  the  Taber  sites  during  August  11-12/76 
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Figure  A-3. 


Potential  evaporation  fluxes,  air  temperatures  and 
wind  speeds  at  the  Taber  sites  during  August  13-14/76 
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WATER  CHARACTERISTIC  CURVES-  SITE  2 
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:igure  B-l . 


Laboratory-determined  water  characteristic  curves  for 
Site  2 
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Figure  B-2.  Comparison  of  laboratory  and  field-determined 
water  characteristic  data  for  Site  2 
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Figure  B-3.  Vapor  equilibration  experiments  using  acid 
solutions  for  Site  2  (0-20  cm) 


APPENDIX  C 


HYDRAULIC  PROPERTIES  OF  POROUS  MATERIALS  FROM  THE  TABER 


SITES 
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graphical 


Appendix  includes  the  following  relationships 
form  as  determined  in  the  laboratory: 
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>cm  H2O 
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Figure  C-l.  ip  (0)  relationship  for  Site  1  (0-30  cm) 
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Figure  C-2.  Sg  (ip)  relationship  for  Site  1  (0-30  cm) 
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Figure  C-3.  Kr  (p)  relationship  for  Site  1  (0-30  cm) 


>  cm  /sec 
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Figure  C-4.  K  (0)  relationship  for  Site  1  (0-30  cm) 
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%  e 
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Figure  C-5.  ^  (0)  relationship  for  Site  1  (90-150  cm) 
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Figure  C-6.  Sg  (ip)  relationship  for  Site  1  (90-150  cm) 
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Figure  C-7.  «r  (ip)  relationship  for  Site  1  (90-  150  cm) 


cm  /sec 
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Figure  C-8.  K  (e)  relationship  for  Site  1  (  90-  150  cm  ) 


cm  HoO 
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Figure  C-9.  ip  (S  )  relationship  for  Site  2  (0-20  cm) 
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Figure  C-10.  Sg  (^)  relationship  for  Site  2  (0-20  cm) 
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Figure  C-12.  Sg  (ip)  relationship  for  Site  2  (150-200  cm) 
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Figure  C-13.  (ip)  relationship  for  Site  2  (150-200  cm) 
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Figure  C-14.  K  (0)  relationship  for  Site  2  (1  50-200  cm) 


APPENDIX  D 


DERIVATION  OF  THE  FINITE  DIFFERENCE  EQUATIONS  OF  THE 


ONE-DIMENSIONAL  COUPLED  NON-I SOTHERMA L  FLOW  USING  A 


CRANK-NICHOLSON  APPROXIMATION  METHOD 


The  set  of  equations  to  be  solved  is  (4-1)  subject  to 
the  initial  and  the  boundary  conditions  (4-2)  and  (4-3) 
respectively.  The  spatial  and  time  domains  are  discretized 
using  the  grid  of  Fig.  4-2.  The  following  modification  of 
the  Crank-N ic holson  method  of  finite  differencing  (Richtmyer 
and  Morton,  1967)  is  used  to  approximate  the  partial 
derivatives  appearing  in  the  above  equations: 


CD-I) 


(D  —  2 ) 


(D-3 ) 


are  arbitrarily  given  by  (4-4r) 


and  (4-4  s)  respecti  vely . 

Using  the  above  approximations  to  equations  (4-1)  they 
become 
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23  b 


f: 

t 


)/At  -fl-  ff1  )/(2Az) 


-k0a  -vL  -W-7  >/<24z> 


.J'!4  .J  .J-'  _J 


+  dTi<-i<4  <T4i  +TC+,  -T; -Tf  ,/<2^> 

-Ti.,  -T>:!  )/(24z)  -<k£^  -4.^1  ] 


(D-4) 


,  j-A 


C(  “  (T^-T2'1  )/4t  =  (1/4z)[3,2.A  Cfi  -t/+,  -Tj'-lf'  )  /  (2Az) 


-J~  Vi.  j  /_•  i  j-i 

(T^T-;  -If.<  -Ti_,  )/{2Az) 

+  DlIvU  <  +4,  +  tin  -  tf  "1  >  /  < : ) 

-< 9^'4  » 3 

j-  »/z 

where  the  terms  of  the  form  k .  are  given  by  (4-4t)  . 
Rearranging  terms  and  making  use  of  ( 4  —  4 f >  ,  the  above 
equations  become 

„  J-Vz  ,  j  ,  0  J-Vz  /T_  .  „  J  “  A  ,  |  J  „  3  “A  j 

■‘tyt-vl  4\-<  t(2Fc  /r+K4w*  Kf.'-'/A  >4r  / -k>p.>j4 4\+, 

+  <474.  T^'  * 

>  +k^4  <  C  -  *{" )  +  -Tf 

-Tf'  )  +  24z<fC&%  -^1*4  )  *  (2F/-‘/i  /r)^' 


(D-5) 


-  gi",/4  Az/2)t{.  -  J 


U,  -°4;-4 
J  .  ^  J'A  J 


(D-6) 

-(4'-4+c/i"!/l'  ^z/2)T(i,  +(2C^‘'/i  /r  +  n4 401,4b  T-f 

r'A  _  > 

l+(/2_ 

+  (Dl"i>V  +D4-4  >4^ 

A'1  >  * 34a (if-!  -if'1) 

♦■>&!*  -i?H  >  +Di;'-ti 

+  c^  (Az/2)  (T^'Z  -T&  )  +(2C/'<A  /t)Tf'  (D-7) 

Using  (4-4c)  to  (4-4q)  ,  the  above  equations  become  (4-4a) 
and  (4- 4b) - 
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Boundar  y  Conditions 


The  top  boundary  condition  for  the  pressure  head 
equation  is 

qw/^  =  ±E  =  -K^  -Dr^T/^Z  +Ki  (D-8) 

which,  after  discretization  becomes 

J'Vz 


si  =  -4’A 

~dTc-*/  (^+Tf'  -t£-,  -T&  )/(24z)  +  4', 4 


(D-9) 

where  i=2  (top  node)  and  the  (i-1/2)  subcripts  are 
approximated  by  i.  Examination  of  (D-4)  indicates  that  the 
underlined  terms  actually  constitute  (D-9) .  Therefore, 
replacing  the  terms  which  are  identical  to  (D-9)  by  term  E 
in  (D-4)  results  in 

j-’4 

^  )  /.ar  -  i  i/azj  l 


rf4  =  (i/az)[k44(4.  +4' -t’-'l'/"  )/<2^> 

+  (4.  ♦*£  >/(2Az)  +Ef  -kC4  ] 


( D- 10) 


which  after  some  rearrangement  becomes 


(2P,^/r*KtJc;4  )4  ^  -Drd^  T‘ 

‘244/r-Ktl>,4)^>'  ♦k4^.V4’  -D-n-Vj  Tc" 

+  2dz  (EC  -  k!"^-  ) 
or,  using  (4-4c)  to  (4-4i) 

-c^  -c^  =  H  (2) 

where  H  (2)  is  the  right-hand-side  of  (D— 11). 


L  +i 

+dj-'a  4: 

rc+a  T^- 

(  D-  1  1 ) 


(D-12) 


If  (D-9)  is  solved  for  the  terms  underlined,  with  i=3 
(assuming  that  E 1  £  -*E )  and  the  result  replaced  in  the 
heat  flow  equation  (D-5),  the  following  result  is  obtained: 
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'"vUt/"'  )/(2Az) 


*  K\  (D-13) 

where  ( i— 1 )  =  2  (top  node).  After  some  rearrangement,  the 


above  equation  becomes 
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n  /K  1j_!4. +_J-K  „  J*' 

+  DLc-*4  /K'f>c-I/J_+C^  *■  ^  ^2/2) 

where  T^_,  =  is  the  specified  surface  temperature,  or 


(  D- 1  4 ) 


using  (4-4c)  to  (4-4p) 

[  B3<  -  (A,' A3 /A,  )  ]T \  -C3'TJ  *(b;-T’)^  -C5'v^  =  H1  (3) 
where  H 1  (3)  is  the  rignt-hand-side  of  (D-14). 


(D-15) 


With  regard  to  the  bottom  boundary  conditions,  the 
water  flux  equation  (2-38)  --  after  discretization  using 

forward  differences  --  becomes 

*KU(vk_+i  /Az  ‘dtMt4I  *Tl)/Az  +Kl  =  ±w  (D-16) 

where  w  is  the  specified  bottom  boundary  water  flux.  Solving 
(D-16)  for  results  in 

'|V+l='fL  -Az(1-W/KL)  -(Dt/K)l  (TLf1  -Tt)  (D- 17) 

while  from  the  bottom  boundary  condition  ( D  —  1 8)  for  the  heat 


flow  equation 
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C3T/Jz)  iz_l  -  (ILt1  -\)/Az  =  0  (D-18) 

one  obtains 


(D-1S) 


APPENDIX  E 


DOCUMENTED  COMPUTER  PROGRAMS 

El  Introduction 

The  present  computer  program  simulates 
one-dimensional  simultaneous  flow  of  water  and  heat  in 
subsurface  (4-1).  It  consists  of  a  main  program  a 
subroutines  or  function  routines.  Most  transport  param 
and  coefficients  are  written  in  separate  subroutin 
order  to  facilitate  changes  by  simply  inserting 
appropriate  program  segment.  The  program  was  initiall 
in  WATFIV  for  debugging  purposes  and  subsequently  in  FO 
G  for  shorter  execution  time.  Double  precision  arith 
was  used  throughout  for  the  various  calculations  involv 
the  model.  The  amount  of  computer  storage  require 
dependent  on  the  size  of  the  finite  difference  scheme  t 
solved.  The  present  program  is  set  up  for  50  nodes.  In 
to  alter  the  size  of  the  system,  the  dimension  state 
must  be  changed. 


The  main  program  is  the  control  program  which  di 
the  sequence  of  operations  for  solving  the  system  of  co 
partial  differential  equations  (4-1).  It  calls 
appropriate  subroutines;  it  controls  the  input  and  o 
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data,  and  boundary  and  initial  conditions;  it  check 
time  increments,  the  number  of  iterations  and  conver 
constraints. 

E2  Subroutines 

Suoroutine  DTCOEF  calculates  the  thermal  diff 
coefficients  D-^  and  Dy  ;  the  non-isothermal  water 
q^/^  I  the  heat  flux  due  to  the  movement  of  water  qcV  an 
total  heat  flux  q^  . 

Subroutine  DPSICF  calculates  the  total 
conductivity  Ko;  . 

Subroutine  HYCOND  adjusts  the  measured  satu 

hydraulic  conductivity  for  temperature  usinq  (4-15).  It 
provides  for  the  calculation  of  the  unsaturated  hydr 
conductivity  function  from  the  non-isothermal  water 
eq  uatio  n  (4- 1 d) . 

Subroutine  EVAPO  calculates  the  evaporation  rate  f 
water  conservation  equation;  that  is,  from  the  time  c 
of  the  water  content  of  the  top  finite  difference 
element  and  the  known  water  flux  at  the  bottom  of 
element.  This  water  flux  is  estimated  from 

non-isothermal  flux  equation  ( 4 -  Id)  employing  the  spec 
hydraulic  conductivity  function  (K  vs  Q  )  for  the 
element . 

Subroutine  HCAPT  calculates  the  volumetric 
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capacity 

of  t  he 

subsurface 

layer  s 

as  functions 

of  their 

mine  ral 

matter. 

organic  mat 

ter  and 

water  fractions 

(3-15). 

Subroutine  HFLBX  calculates  the  water  vapor 
conductivity  ;  the  conductive  gcj  and  latent  g^  heat 
fluxes;  the  vapor  thermal  diffusion  coefficient  D-^  and  the 
vapor  flux  gv/^. 

Subroutine  FCOEF  calculates  the  specific  water  capacity 
(2-47)  . 


Subroutine  DLCOEF  calculates 
for  heat  transfer  Dj_  (4-1c). 

Subroutine  TRIDAG  (Carnaha 
system  of  linear  simultaneous  egu 
matrix;  such  a  system  results 
difference  scheme  used  here.  The 
IF  through  L  and  their  su 
superdiagonal  coefficients  are 
and  CZ.  The  computed  solution  vec 
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9  6  9)  sol  ves  a 
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plicit  finite 
numbered  from 
iagonal  and 
arrays  AZ,  3Z 
in  the  array 


VZ. 


Subroutines  SPLINE  and  CALCCF  and  function  PCDBIC  form 
the  piecewise-cubic  interpolation  program  (Conte  and  de 
Boor,  1972).  This  cubic  spline  program  is  called  extensively 
to  interpolate  input  data  for  a  number  of  parameters 
reguired  in  the  program. Subroutine  SPLINE  uses  Gauss 
elimination  adapted  to  take  advantage  of  the  tridiagonal 
character  of  the  resulting  coefficient  matrix  (Conte  and  de 
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Boor,  1  972)  to  calculate  the  spline  coefficients  C(2,i)  ,  i  = 
2r  /  N,  given  the  numbers  C(1,i)  =  f(Xj_),  i  =  1,  ...  , 

N+ 1 ,  C(  2,  1)  =  f*  (X|_)  and  C(2,  N  +  1)  =  f  ’  (x^,  )  ,  where  f  ( x  ^ ) 
comprises  the  input  data  corresponding  to  x^  and  C(2,1)  and 
C (  2 ,  N+ 1 )  are  reasonable  approximations  to  the  end  slopes  of 
f  (XjJ  .  Subroutine  CALCCF  calculates  the  spline  coefficients 
C(3,i)  and  C(4,i),  where  i  =  1,  ...  ,  N.  Function  P  CUBIC 

evaluates  the  cubic  spline  equation  (E-1) 

f  (x-  )-.C  (1,i)+C  (2 ,  i)  (x-xj  +C  (3  ,  i)  (x-xc)  2+c  (4,i)  (x-xc)  3  (E-1) 

for  any  particular  point  once  the  spline  coefficients 
C(j,i),  j  =  1/..-4,  i  =  1,...N  are  known. 

Subroutine  KSPLEN  interpolates  the  specified  hydraulic 
conductivity  function  (K  vs  0 )  data  using  cubic  splines  and 
adjusts  the  conductivity  values  for  temperature. 

Subroutines  DENS,  VISCO,  VAPSPL,  CLQCF  and  DA  If  C  F 
calculate  water  density  ^  ,  viscosity  y  ,  latent  heat  of 
vaporization  of  water  L',  saturated  vapor  pressure  pj  , 
specific  heat  of  liquid  water  c ^  ,  and  molecular  diffusion 

coefficient  of  water  vapor  in  air  Da  ,  respectively,  as 
functions  of  temperature.  All  the  above  subroutines  call  tne 
cubic  spline  interpolation  program  to  interpolate  the  input 
da  ta . 


Subroutines  THESP2  and  THCOND  calculate  the  volumetric 
water  content  and  thermal  conductivity  as  functions  of  the 
negative  pressure  head  and  water  content,  respectively.  Both 
subroutines  interpolate  the  input  data  using  the  cubic 
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spline  interpolation  program  mentioned  above. 
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E3  List  of  Variables* 


cc ,  cco 

volumetric  heat  capacity  of  porous 

material  (C)  at  current  and  previous 

time  steps 

CL Q,  CL  QO 

specific  heat  of  liquid  water  (c^) 

at  current  and  previous  time  steps 

CSATi,  i=1 , ..  .  ,  4 

saturated  hydraulic  conductivity 

(Ks«t  )  of  layer  i 

DAW 

molecular  diffusion  coefficient  of 

water  vapor  in  air  (Da) 

DELT 

time  increment  (At) 

DELZ 

depth  increment  (Az) 

DP  SI ,  DPSIO 

total  water  conductivity  (K^)  of 

current  and  previous  time  steps 

DP  SI V,  DPSIVO 

vapor  conductivity  (Kvpv)  of  current 

and  previous  time  steps 

DT,  DTO 

thermal  water  diffusivity  (D-j-)  of 

current  and  previous  time  steps 

DTHETA 

water  diffusivity  (Dg) 

DTH  V 

vapor  diffusivity  (Dy  ) 

DTL 

liquid  thermal  diffusivity  (D -^  ) 

DT V ,  DT  VO 

thermal  vapor  diffusivity  ( D-y. )  of 

current,  and  previous  time  steps 

EPS 

a  small  positive  number  (e) 

FF ,  FO 

specific  water  capacity  (2>0/cAp)  or 

1  Units :  [length]  =  cm;  [time]  =  hr;  [mass]  =  gm ; 

[temperature]  =  C;  [heat]  =  cal. 
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GR 

HVAP 


specific  storage  (Ss ) 

acceleration  of  gravity  (g) 

latent  heat  of  vaporization  of  water 

(L') 


J 

time  step 

index 

K,  KO 

hydraulic 

conductivity  (K) 

current  a 

nd  previous  time  steps 

KM  AX 

maximum 

number  of  iterations 

time  step 

KODE  1 

specifies 

the  top  boundary  cond 

of  evapor 

ation  or  infiltration 

K0DE2 

specifies 

after  how  many  time 

the  outpu 

t  is  to  be  printed 

K0DE4 

specifies 

if  the  coupled 

uncoupled 

equations  are  to  be  u 

K0BE5 

specifies 

the  bottom  bou 

condition 

of  recharge  or  discha 

KOUNT 

count  ind 

ex 

KS  AT 

saturated 

hydraulic  conduct 

(Ksat  ) 

L 

total  num 

ber  of  space  nodes 

LAM  DA,  LAMDAO 

apparent 

thermal  conductivity  (J 

current  a 

nd  previous  time  steps 

MI 

dynamic  v 

iscosity  of  water  (ja) 

MIT REF 

dynamic 

viscosity  of  wate 

reference 

temperature  {j4(T0)} 

MO 

node  belo 

w  which  the  mineral 

or 


MW 

Ni 

N 1 

N2 

N3 

N5 

N7 

N8 

N9 

NI  1 

N6  6 

NJ 

NR,  NI,  Nil 

NS  TOP 

Pi ,  i= 1,2, 3, 5, 7, 
8,9,11,66 
Pi  (1,1)  ,  1=1,  Ni 

PSIO,  PSI,  PSII 

PSI66 


organic  mat 
constant 
molecular  wei 
number  of  spe 
various  param 
number  of  ^ 

n  ii  y 

it  ii  l  i 

"  "  K 

"  »  p5 

rv 

"  "  Ct 

"  "  D* 

"  "  3 

ii  it  ^ 

number  of  tim 

depths  at 
hydraulic  con 
upper  limit  o 
increments 
end  slope  of 
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ter  fractions  remain 

ght  of  water  vapor  (M) 
cified  input  data  of 
eters 

(T)  input  data 
(T)  "  " 

(T)  "  " 

(0> 

(T)  "  " 

(T )  "  " 

(T)  "  " 

(0 ) 

(0 )  "  " 

e  increments 
which  the  saturated 
ductivitv  changes 
n  the  number  of  time 

input  data  f  (x^ ) 


input  data  corresponding  to  Ni  (see 


Ni) 


pressure  head 
current  and  next 
input  pressure 


«j») 

time 

head 


of  previous, 
steps 

data  in  the 


water  characteristic  function 
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QH 

total  heat  flux  (g^) 

QCD 

conductive  heat  flux  (qcc/) 

QCV 

convective  heat  flux  (gcv ) 

QS 

latent  heat  flux  (q^) 

QW,  QWO 

water  flux  (dv/^)  of  current 

previous  time  steps 

and 

QVAP 

vapor  flux  (qv/^  ) 

R 

universal  gas  constant  (R) 

RO 

liquid  water  density  (^) 

ROTREF 

liguid  water  density  at 

reference  temperature  {  (>(T0  )  } 

the 

SL  total  length  of  subsurface  profile 

studied  (L) 

SR H  soil  relative  humidity  (h) 

SS  specific  storage  (S5  ) 

SVP  saturated  vapor  pressure  (p*) 

TO,  T,  TI  porous  medium  temperature  (T )  of 

previous  current  and  next  time  steps 
Ti (I) ,  1=  1 , Ni  input  temperature  data  corresponding 

to  Ni  (see  Ni) 

THETAO,  THETA,  THETAI  moisture  content  (0)  of  previous, 

current  and  next  time  steps 

THETA5  input  moisture  content  data  in  the 

M0)  relationship 

THET11  input  moisture  content  data  in  the 

rela  tionship 


TINCR 


incremental  time  step 
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TREF 

reference  temperature  (Te  ) 

TS 

ground  surface  temperature 

W 

specified  water  flux  at  bottom 

boundary 

XM 

volumetric  mineral  matter  fraction 

XMAXST 

maximum  simulation  time 

XO 

volumetric  organic  matter  fraction 

(x0) 

Z 

incremental  depth  (z) . 
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E4  Program  LI  sting 


1 

2 

3 

4 

5 

6 
7 
3 
9 

to 
1 1 
12 

13 

14 

15 

16 
17 
16 

19 

20 
21 
22 

23 

24 

25 

26 
27 
23 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 
4  3 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 


IMPLICIT  REAL*8 (A-H,0-Z) 

COMHON/D5WHE/  PSI  (50)  , THETA (50)  ,T(50)  ,QW(50)  ,QH(50) .DELTAT (50)  , 
1LAMDA (50)  ,CLQ (50) ,SV  P(50)  ,30 (50) ,MI (50)  ,P2  (4,50) ,T2(50) , 

2T1  (50)  ,  PI  (4,50)  ,T3(50)  ,P3(4,50)  ,  Tit  ST  11  (50)  ,P11  (4,50)  ,H7AP(5Q)  , 
3TH£TA5(50)  ,P5(4,5(J)  , T9  (50)  ,P9(4,50)  ,FF(50)  , 

4PSI6  6  (50)  ,P6  6  (4, 50)  , T7  (50) ,P7(4 ,50)  ,T8  (SO)  , 

5P8 (4,50)  ,  SRH  (50)  ,K (50)  ,QWO(50)  ,DAU (50)  ,  OT  (50) , 

6THETA0 (50) ,NJ, L, LI  , K L, 3W , DELZ 
C 

COHMON/NN/  N1,N2,N3,N5,N7,N8,N9,M11,N66 

DIMENSION  PSIO  (50) ,TO (50) ,  LAMD A 0 (5 0) , APS (50) , DTL ( 50)  , 

1CP8 (50) , BPS  (50)  , CC (50)  ,CCO  (50)  , D L ( 50 )  , DLO  ( 50) , QV AP  (50) , 

2CBARPR  (50)  , BBAHPR  (50)  ,DTO(50)  , A BAR  (50)  ,CBAR (50)  , 

3 A (50) ,C  (50)  , B (50)  , DP SI  (5  0)  , DPS 10  (5  0)  ,K PLUS (50)  ,DTHETA (50) , 

4KO (50)  ,KSAT  (50)  , DTV  (50)  ,DTVO(50)  , DTH V (50)  , DPS IV (50) , 

5D (50, 3)  ,  E (50, 3)  ,  F (50 , 3)  ,G (50 , 3) , CTT  (SO)  . FTT (SO)  , 

6H (50) ,HPR (50) ,  0  (50) , V  (50) ,TS  (6  00) , £S (50) , II (50) , 

7DELTRF (50)  , KMI N U S  ( 50 )  , TI  ( 50)  ,PSII  (50)  , IKETAI (50)  , 

8XM  (50)  ,  XO  (50)  , ABA  RPR  (5  0)  , 3 BAR (50)  , QCD  (50) ,P  (4, 50) , 

9DPSIVO  (50)  ,FO  (50)  ,CQ  (50)  ,Q»3  (50)  ,CLQO (50)  , 

>DI  (50)  ,  DII  (50)  , Dill  (50)  ,01  (50)  ,GII  (50)  , GUI  (50) 

C 

REAL *8  K ,MW ,  LAMD  A, KS AT  ,  MI , K ITRE? ,H V AP, KMI N US,  KP LOS , LAM DAO  ,  SO 

C - BEAD  1 N  AND  PRINT  PHYSICAL  PARAMETERS  FOR  THE  SYSTEM 

READ  (5, 501)  R , 3 W , SS , SL , GR ,T REF , ROT R  EF , MIT S E? 

WRITE (6, 200)  R,HW,SS,SL,GR,TREF,RUTREF,RITSEF 
200  FORMAT  £ 1H1 ,' LIST  OF  PH  YSICAL  PARAMETERS  POR  THE  SYSTEM*,// 

1 1H0, *GAS  CONSTANT  (  ERG/D EG- C- MO L)  * 

2  ,F10.  1 ,  /  I  HO , ' MOLECULAR  WT.  OF  WATER  (  G3/MOL  )*,P10.1,/ 

31 HO, 'SPECIFIC  STORAGE (CM**- 1) • ,5X,S15.  7,/ 

41H0, • SECTION  L  EN  GTH  (  CM  )  • , 1 4X , ? 1 0- 1 , / 

61H0,  'ACCELERATION  OF  G RA VIT Y (C M/SEC**2 )  » , F 10. 3 ,/ 

71H0,  'REFERENCE  TEMP.  (DEG.C)',12X,F10.3»/ 

81H0, 'WATER  DENSITY  AT  REF . TEMP.  (GM/CC)  • , F  1  1. 5 . / 

91H0, ' WATER  VISCOSITY  AT  REF.  TEMP,  (GM/CM-HR)  F6. 3,////) 

501  FORMAT (2  P 1 0. 1, DIO. 5,  FI  0-  1,2P1Q.3,/2?10.5) 

C - READ  IN  AND  PRINT  NUMERICAL  INFORMATION 

READ (5, 502) NSTOP,L,NJ, XMAXST,EPS, KM AX, DELZ 
50  2  FORM  AT (316, F3-2,F7. 5,14,510.3) 

J  N  J=  N J 
L 1=L- 1 

DELT  =  XM A XST/XN  J 
F ATIO=  DELT/DEL  Z 
RATIO 3  =  RAT IO/D  ELZ 

WRITS  (6,  2  30)  XMAXST,DELT,NJ,L, DELZ, EPS, KMAX 
230  FORMAT (1H1,' NUMERICAL  IN  FOR M ATION • , //I  HO , • M A XI  MOM  SIMULATION  TIME 
1  (HR)  ',  1X,G15.7,/1K0, 'SIZE  OF  TIME  INCR EMENT (HR )  • , 3 X , G 1 5- 7, /I HO , 

2 '  N  UMBER  OF  TIME  INCR  EM  ENTS  •  ,1 1  5  ,/1  HO  , 

3'NUMBER  OF  SPACE  NODES'  ,ri9,/1H0, 

4  *  SIZE  OF  SPACE  INCREMENTS {  CM  )',  F 1 0. 3 , /I  HO , • EPS ILON * , 23X , F 10. 4, 
5/1H0, ' MAXIMUM  NUMBER  OF  ITERATIONS ', 2X , 110 ,////) 

C - READ  IN  AND  PRINT  BOUNDARY  CONDITIONS 

TINCR=0. 0D0 
DO  81  1=1, NJ 

TS (I) =24. 4D0*7. 7D0*DSIN  ( 6. 2 83D0 *TI NCR/24- 0DQ*O_ 392 7D 0) 

81  TINCR=TINCR+DELT 

IF  (NJ- LE-  1000)  GO  TO  93 

C -  PICEUICE  LINEAR  SURFACE  TEMPERATURE  DISTRIBUTION 

TS (1) =32- 5DO 
DO  82  1=2,7 


61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

30 

81 

32 

33 

84 

35 

36 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

1  10 

111 

112 

113 

114 

115 

1 16 

117 

118 

119 

120 
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82  TS  (I) =TS  (I-  1) *0.  13666D0 
DO  83  1=8,14 

83  TS  (I)  =TS  ( X  —  1 )  *0.  05571D0 
DO  84  1=15,21 

84  TS  (I)  =T3  (I-  1) -0.  53D0 
DO  85  1=22,28 

85  TS (I) =TS  (I-  1) -0.  59286D0 
DO  86  1=  29,35 

86  TS (I) =TS (I-  1) -0.  4957  IDO 
DO  88  1=36,42 

88  TS  (I)  =TS  (I-  1) -0.  13143D0 
DO  89  1=43,56 

89  TS(I)  =TS  (I-  1)-0.  19692D0 
DO  90  1=57,63 

90  TS(I)=TS  (I-  1)-0.  13333D0 
DO  91  1=64,70 

91  TS(I) =TS  (I- 1) +0.  22833D0 
DO  92  1=71,77 

92  TS(I)  =TS  (1-  1)  *0.  755D0 

93  READ  (5, 121)  KODE 1 , KO DE2 , KODE3 , KODE 5 

121  FORMAT  (413) 

C- - BOTTOM  BOUNDARY  HATER  FLUX 

IF  (h'ODE  1 )  5,5,10 

5  QH  (L) =0 - 0D0 
GO  TO  122 
10  QH  (L) =  0.07D0 

122  IJK=NJ/2 

DO  94  J= 1 , I JK 
JJ=J*I JK 

94  HBITE{6, 141)  J ,TS  (J)  ,TS ( JJ)  ,  JJ 

141  FORMAT ( 1  H  , 13 , 3X , • SO RFAC E  TEM PE B AT URE ( DEG. C)  • , G  15. 7, 3X , G 1 5.  7, 

1  2  X  ,  1 3 ) 

S=QW  (L) 

C — — READ  IN  AND  PRINT  OUT  INITIAL  CONDITIONS 
READ  (5,503)  (PSI(I)  ,1=2, L) 

READ(5,503)  (T(I),I  =  2,L) 

503  FORMAT  (  6F10. 3) 

T  (  2)  =TS  (  1) 

241  8 RITE  (6 , 240 ) 

240  FORMAT ( 1  HI ,’ LISTING  OF  INITIAL  CONDITI ONS ’ , // 1  HO , 1  OX , * INCREM ENT  *  , 
1  18X,  ’POSITION'  , 17X, • PS  I* ,  5X, 'MOISTURE  CONTENT (THETA)  ' ,5X, 

2 ' TEMPE RATUR E ' ,/1H  , 4  IX  ,  •  (CM )  • , 1 2X , '  (CM  U20)  • , 1 8X , •  (DIM. )  • , 

313X, ’  (DEG. C)  ' ,///) 

C 

KL=L* 1 

READ  (5, 281)  N  1 ,  N  2,  N  3  ,  N  5,  N7,  N8,  N  9 ,  N  1 1 

281  FORMAT  (815) 

READ  (5, 2222)  (T 1  ( I)  ,  P 1  ( 1 , 1)  ,  1=  1  ,  N  1 )  ,  P  1  (2,  1)  ,P1  (2,N1) 

READ  (5, 2222)  ( TH  ET1  1  (I )  ,  P 1 1  ( 1 , 1)  ,1=  1 ,  N1 1)  ,  P1 1  ( 2 ,  1)  ,P  1 1  (2 ,  Hi  1) 

READ (5,282)  N66,MO 

282  FO RM AT  (215) 

C- —  CHARACTERISTIC  CURVES 

READ  (5, 3333)  (PSI66(I},P66(1,I)  ,  1= 1 , N6 6 )  , P6 6  (2,  1)  , P 66 ( 2 , N 66 ) 

3333  FORMAT  (D15. 7, r  10.4) 

READ (5,154)  NR  ,NI,NII,CSAT1 , CS AT2 , CSAT3 ,CSAT4 
154  FORMAT  (315,  4F10.  5) 

DO  1  1=2, NR 
1  KSAT (I ) =CS AT  1 
NN 1=  NR ♦  1 
DO  2  1=  N  N  1  ,  NI 
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136 
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139 
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141 
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143 

144 

145 
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148 

149 

150 
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152 
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154 
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2  KS  AT (X) =CS AT2 
NN  2=  NI ♦  1 

DO  3  I=NN2,L1 

3  KSAT(I) =CSAT3 

4  KSAT (L) =  CSAT  4 

CALL  THESP2  (P 51 , K L, TH ETA , N 66 , PSI66 , P6 6 , L) 

C - 

DO  6  1  =  2, L 
IF(I.GT.MO)  GO  TO  7 
REAO  (5, 1  111)  XM  (I)  ,XO  (I) 

GO  TO  6 
7  XM (I) =  XM (MO) 

XO  (I)=XO  (MO) 

6  CONTIS  UE 

1111  FORMAT (F10. 3,F10. 5) 

HEAD (5,2222)  (T9  (I)  ,  P9  ( 1  , 1)  ,1= 1 , N9)  ,  P9  (2,  1 )  , P9  (2,  N9) 

HEAD (5, 2  2  22)  (T3  (X)  , P8 ( 1 , 1)  , I=1,N3) ,P3  (2,  1) ,P8  (2,N8) 

READ  (5,22  22)  (T7  (I)  ,P7  (1,1)  ,1=1,  N7)  ,P7  (2,  1)  ,P7  (2,N7) 

READ  (5, 2222)  (T3(I),P3(1,I)  ,I=1,N3),P3(2,1),P3(2,N3) 

2222  FORMAT  (2F10.  5) 

READ (5, 2  2  22)  (T2  (I)  , P2  ( 1 , 1)  , 1= 1 , N2 )  ,P2  (2,  1 ) , P2  ( 2, N 2) 

READ (5,2  225)  (THETA 5  (I)  , P5  (1, 1)  , 1= 1 , N5) , P 5  (2 , 1 )  ,P5 (2,N5) 
2225  FORMAT  (F 10. 5, DIO- 6) 

DO  250  1-2, L 
2  =  1-2 
Z=DELZ*Z 

211  KBIT E(6, 260)  I , Z , PS  I  (I )  , THET A ( I )  , T (I) 

250  CONTINUE 

260  FORMAT  (120, 6 X , ?2 0 . 8 , 4X  ,3F20 . 3) 

T(KL)=T(L) 

PSI (KL) =  PSI  (L)  -0  2LZ* ( 1 . 0  DO*  9/KS  AT (L)  ) 

KODE4=KO  DE  2 

C - BEGIN  TIME  ITERATIONS 

40  DO  270  J= 1 , N J 

IF  (NJ.GT.NSTOP)  GO  TO  999 
KOUNT=0 

IF  (J.EQ. 1)  GO  TO  693  1 
GO  TO  699 
6981  DO  697  1=2, L 

THETAO  (I ) =TH ET A  (I) 

TO (I) =T  (I) 

PSIO  (I)  =PSI  (I) 

697  CONTINUE 

C - CALCULATE  TRAN5POBT  COEFFICIENTS 

699  CALL  dCA PT (XM , XO , L, THETA ,CC ,KL) 

IF  (J.GE. 43)  GO  TO  999 

CALL  THESP2  (P SI , K L, TH ET A , N6 6 , PSI 6 6 ,P6 6 , L) 

DO  5325  1=2, L 
IF(J.EQ.I)  CCO(I)=CC(I) 

DELTAT  (I) =DABS  (T  ( I ♦ 1 ) —  T  (I) ) 

CTT  (I)  =  (CC  (I) +CCO  (I)  J/2.0D0 
DELTEF (I) =  D A  3S  (T  (I) -TH  EF) 

IF  (KODE3- EQ.  0)  DL(I)=0.0D0 

IF (KODE3. EQ. 0)  CLQ(I)=O.ODO 
359  IF  (PSI (I) .GE. 0. 0D0)  GO  TO  361 
EN  =  M»'*GS*PSI  (I) 

ED=  R *  (T  (I)  *-2  73.  16D0) 

RE=EN/ED 

SRH  (I) =DEXP  (RE) 

GO  TO  5325 
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216 
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361  SRH  (I) =  1 . 0  DO 
5325  CONTINUE 

698  CALL  VAPSPL  (T ,  L ,  H V A P , N 3 , T 3 , P3 , KL) 

CALL  DENS  (T, L, RO , N 1,T1 , P 1 , KL) 

CALL  THCOND(L, LAMDA,  PSI  ,  TH  ET  A  ,  K  L  ,  XH  ET 1  1,P1  1,N1  1) 

CALL  7 1  SCO  (T,L,MI,N2,T2,P2,KL) 

CALL  SVPSPL (T , L , S  VP  , N7 ,  T7 , P7 , KL) 

CALL  DAWCF  (T,L  ,D  A W  ,  N 9 ,  T9  , P9  ,  KL) 

CALL  FCQEF  (L, PSI, THETA, SS,FF,KL) 

IF  (KODE3. EQ. 0)  GO  TO  605 
CALL  CLQCP  (T,L ,CLQ,N 8, T8, P8, KL) 

605  CALL  K3PLEN (THETA, K, N5 ,THETA5, P5,L , PSI ,T, MI ,RO, 

1T1 ,T2, PI ,P2, N1 ,N2,KL,KSAT) 

CALL  HFLUX  (  DTH  V , DT V , Q V AP , DPS I V , QS ,  QCD, KODE3) 

CALL  EVAPO  (EV,  DE LT  ,  DTV  ,  D PSI  V ,  DT  L,  J  ,  KOU  NT,  KODE 3  ,  KOD  E5) 

606  CALL  HXCOND ( J , KS AT , MIT  REF , ROTR EF , DTH V , DPSI V , 

1DELTRF,DTV, EV) 

CALL  DPSICF  (DT HV , DP  SI , KS AT, DTH ETA , D PSI V) 

CALL  DTCOEF (DTV,DTL,DPSI,EV,NR,NI,QCD,QS,KODE3) 

IF  (KODE3. EQ.O)  GO  TO  696 
CALL  DLCOEF (DTHV ,DL,DPSIV) 

696  IF  (KOUNT.EQ. 0)  GO  TO  777 
778  KOUNT=KOUNT+ 1 

IF  (KODE2. EQ. 1)  GO  TO  6011 
IF  (J-K0DE4. N E. 0)  GO  TO  341 
6011  WRITE (6,6768) 

6768  FO  RMAT  (1  HO,  10  X  ,  •  THET  A*  ,  1  4X,  '  K*  ,  12X,  •  QH  *,11X/'CC  MSX.'QV, 
112X, ' DTL • , 1 2X ,  'PSI* , 12 X,  'T* ) 

DO  6769  1=2, L 

6769  WRITE  (6, 6770)  I,  THETA  (I )  ,  K  (I)  ,  QH  (I)  ,  CC  (I)  ,  QH  (I)  ,  DTL  (I)  ,  PSI  (I)  , 

IT  (I) 

6770  FORMAT  (1  H  ,  I  2,  21 , 8G  1  5.  6) 

WRITE (6, 6767) 

6767  FORMAT  (1H0,  10X,' THETA*  ,  12X, 'DT' , 12X, 'DTV*  ,  12X, '  DL ' ,10X, 

1  *  L  AMD A  '  , 1 1 X  ,  ' DP  SI  V ’  ,1 IX, 'DPS I' ,1QX,'QVAP*) 

DO  5163  1=2,30 

5163  WRITE (6, 5164)1 , THETA  (I)  , DT (I)  ,DTV  (I)  ,DL(I)  , LAMDA (I) , DPSIV (I)  , 

1 DPSI (I)  , QVAP  (I) 

5164  FORMAT  (1  H  ,  I  2 , 2X  ,  8G  1  5.  6) 

341  DO  340  1=2, L 

FTT  (I)  =  (FF(I)+FO  (I)  )  /2-ODO 
CQ  (I) =  (CLQ  (I) +  CLQO  (I) ) /2.QDQ 
QWW  (I)  =  (QW  (I)  +QWO  (I)  ) /2-0D0 
340  CONTINUE 

DO  720  1=3, LI 

APR(I)  =(LAMDA(I)  +LAMDA  (I  —  1 )  +LAMDAO  (1-1) *  L  AMD AO  (I) ) /4.QD0 
1  +CQ  (I)  *RO  (I)  *Q W W  (I)  *  DE  LZ/2.  0D0 
718  CPR  (I)  =  (LAMDAO  (I)  +  LAMDAO  (1+  1)  «■  L  AMD  A  (1+1)  *LAMDA  (I)  )/4.0D0 
1  — CQ (I)  *RO (I) *Q WW  (I)  *  DEL2/2. ODO 
716  3PR  (I)  =  (  (CC  (I)  «-CCO  (I)  )  /RATIOR)  *-APR  (I)  ♦  CPR  (I) 

720  CONTINUE 

APR  (2)  =  (LAMDA  (2)  LAM  DAO  (2)  )  /2.  ODO  +  CQ  (2)  *RO  (2)  *QWW  (  2)  *DE  LZ/2.  ODO 
CPR (2)  =  (LAM  DAO  (2)  +L AMD  A  (2) >L AMD  AO  (3) *L AMD A  (3) ) /4. ODO 
1-CQ  (2)  *RO  (2)  *QWW  (2)  *DE LZ/2.  ODO 
BPR  (2)  =  (  (CC  (2)  +CCO  (2)  )  /RATIOR)  ♦  A  PR  {  2)  +CPR  ( 2) 

APR  (L)  = (LAMDA (L)  ♦ LAM DA  (LI) ♦ LAM DAO (LI)  +  LAM DAO (L)  )  /4 . 0  DO 
1 *CQ  (L)  *3  0  (L)  *Q WW  (L)  *  DE LZ/2. ODO 

CPR  (L)  =  (LAMDA  (L)  *LAMDAO  (L)  )  /2.  ODO-CQ  (L)  *80  (L)  *QWW  (L)  *DELZ/2.0D0 
BPR  (L)  =  (  (CC  (L)  *CCO  (L)  )  /RATIOS)  +  APR  (L)  *CP3  (L) 

DO  740  1=3, LI 
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ABARPR  (I )  =  (DL  (X)  +  DL  (I-  1)  +  DLO  (1-1)  ♦DLO(I)  )  /4.QD0 
738  CBAHPR  (I)  =  (DLO  (I)  +DLO(I+1) ♦  D  L  (I  ♦ 1)  +  DL ( I) ) /4. 0  DO 
736  BBABPR  (I) =ABARP3 (I)  +CBARPR (I) 

740  CONTINUE 

ABARPB  (2)  =  (DLO  (2) +  DL  (2) ) /2. 0D0 

CBARPR  (2)  =  (DLO  (2)  ♦  DL  (2)  +  DLO  (3)  +  D  L  (  3 )  ) /4. 0  DO 

BBARPB  (2) =ABARPB  (2)  +C3ABPa(2) 

A3ARPR  (L)  =  (DL  (L)  +DL (LI) +DLO (LI)  +DLO  (L)  J/4.QD0 
CBARPB  (L) =  (DLO (L) +DL (L) ) /2. ODO 
BBABPR  (L) =ABAR?H  (L)  +CBARPB (L) 

732  DO  760  1-3, LI 

ABAS  (I)  =  (DT  (I)  +DT  (I-  1)  + DTO (1-1)  +  DTO  (I)  ) /4 . 0 DO 
758  C B AR (I ) =  (DTO  (I)  +  DTO  (1+1)  +  DT  (1+ 1 ) +DT  (I) ) /4.0D0 
756  BBAB  (I) = ABAR  (I) +CBAB  (I) 

760  CONTINUE 

ABAB  (2) =  (DT  (2) +  DTO (2) ) /2.0D0 

CBAR  (2)  =  (DTO(2)  +  DT ( 2 )  +  DTO  (3)  +DT  (3)  )  /4.  ODO 

BBAR  (2)=ABAR  (2) +  CB  A  R  (2) 

ABAR  (L)  =  (DT  (L)  +  DT (L 1 ) + DTO (L 1 ) + DTO  (L) ) /4. 0 DO 
CB AS  (L ) =  (DTO  (L)  +  DT (L) ) /2.0D0 
BBAR(L)=ABAR  (L)  +  CBAR  (L) 

DO  780  1=3, LI 

A  (I)  =  (DPSI  (I) +  DP  SI  (1-1)  +DPSIO(I-1)  +  DPSIO (I ) ) /4.0D0 
768  C (I)  =  ( DPSIO  (I)  +DPSIO  (1  +  1)  +  DPSI  (1  +  1) +  DPSI (I) ) /4.0D0 

76  6  B  (I)=  (  (FF  (I)  +FO  (I)  ) /RATIOR)  +  A (I) +C  (I) 

780  CONTINUE 

A  (2)  =  (DPSI  (2)  +  DPS  10  (2)  )/2.0D0 

C  (2)  =  (DPSIO  (2)  +  DPSI  (  2)  +  DPS  10  (3)  +  DPSI  (3 ) ) /4. ODO 
B  (2)  =  (  (FF  (2)  +  FO  (2)  )  /RATIOR)  +  A  (2)  +  C  (2) 

A  (L)  =  (DPSI  (L)  +  D  P  S I  (LI)  +  DPSIO(L1)+DPSIO(L))/4,ODO 
C  (L)  =  (DPSIO  (L)  +  DPSI(L))/2.0D0 
B (L)  =  ( (FF  (L)  +FO (L) ) /RATIOR)  + A  (L)  *C  (L) 

DO  800  1=3, LI 
IP(I.GE.NR)  GO  TO  775 

KMINUS (I)=  (K  (I)  +  K  (I- 1)  +  KO (1-1)  +  KO (I) ) /4.0D0 
7  88  KPLUS  (I) =  (KO  (I)  +KO  (1  +  1 )  +  K (I  +  1 ) + K (I ) ) /4 . 0 DO 
GO  TO  800 

77  5  KMINUS  (I) =4. ODO/ (1. ODO/K  (I)  +1. ODO/K  (1-1)  +  1  -  0D0 /KO  ( I- 1)  ♦  1-0D0/ 

1 KO (I) ) 

KPLUS  (I)  =4. 0D0/(  1  .ODO/KO  (I)  ♦  1.  ODO/KO  (I*  1)  ♦  1.0D0/K(I+  1)  +1.0D0/ 
1K(I)) 

800  CONTINUE 

KMINUS  (2)=  (K  (2)  +KO  (2)  )  /2.0D0 

KPLUS  (2)  =  (KO  (2)  +  K(2)  +  KO(3)  +K(3)  ) /4 . 0  DO 

KMINUS (L)  =4.0 DO/ (1. ODO/K  (L)  +1. ODO/K (L 1)  +1. ODO/KO (L 1)  +1.0D0/ 

1  KO (L) ) 

KPLUS (L)  =  2. 0D0/ ( 1 .ODO/K (L)  ♦ 1 . OD 0 /KO (L)  ) 

C - FORMULATE  SUBMATRICES  D,  E,  F,  G 

811  D (3, 1) =0. ODO 

D  (3,2) =BPR  (3) -ABARPR  (3)  *  ABA  R  (3) /A  (3) 

D  (3,3)  =-CPR  (3) 

DO  300  1=4, LI 
D (I,  1)  =- APR  (I) 

D  (I  ,  2)  =BPR  (I) 

D  (1,3)  =-CPK  (I) 

300  CONTINUE 

D  (L,  1)  =- APR  (L) 

D  ( L  ,  2)  =B  PR  (L) 

D (L, 3)  =0 . 0  DO 
G (2, 1)  =0. ODO 
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G(2,2)=B  (2)  -  A  (  2) 

G  (2 , 3)  =-C  (2) 

DO  330  1  =  3, LI 
G (I, 1)-~ A  (I) 

G  (1,2)  =B  (X) 

G  (1,3)  =-C  (I) 

330  CONTINUE 

G  (L,  1 )  =  -  A  ( L) 

G  (L  ,  2)  =8  (L) 

G  (L,  3)  =0-  0  DO 
E  (3, 1)  =0. 0D0 
E  (3,2) =C  3 A  RPR (3) 

E ( 3 , 3) =  -C 3A RPH  (3) 

DO  310  1=4, LI 
E  (X,  1)  =-  A  BAB  PR  (I) 

E(I,2)=B3ARP8(I) 

E  (1,3)  =-CB AHPB  (I) 

310  CONTINUE 

E (L, 1)  =- ABARPR  (L) 

E  (L,  2)  =  B8ABPH  (L) 

E  (L,  3)  =0. 0D0 
F (2,  1)  =0 . 0  DO 
F  (2,2)  =C B A B  (2) 

F  (2, 3) =-CBAB (2) 

DO  320  1=3, LI 
F  (I,  1)  =- ABAR  (I) 

F  (I  ,  2)  =BBAR  (I) 

F  (1,3)  =-CBAR  (I) 

320  CONTINUE 

F  (L,  1)  =-A3AB  (L) 

F  (L,  2)  =B 3 A R  (L) 

F  (L,3)=0.GD0 

C - APPLY  TOP  B.C. 

TI  (2 )  =TS  (J) 

C - APPLY  BOTTOM  B.C. 

322  TI (KL) =T  (L) 

30  PSI  (KL)  =  PSI  (L)  -DELL*  (1 . 0D0+ W/K  ( L)  )  -DT  ( L)  *  (T  (KL )  -T  (  L)  )  /K  (L) 

PSII  (KL)  =?SI  (KL) 

IF  (PSI  (L)  .  LT.  0-0  DO)  GO  TO  887 
C 

35  DO  350  1=3, L 

H  (I)  =-  (C  (I)  *  (PSI  (I)  -  PSI  (1  +  1)  )  ♦  A  (I)  *  (PSI  (I)  -PSI  (1-1)) 

1  +C  B A  R  (I)  *  (T  (I)  -T  (1  +  1)  )  +  ABAR  (I)  *  (T  (I)  -T  (1-1)  ) 

2-  (KHINUS  (I)  -KP  LUS  (I)  )  *2.0D0*DELZ-  (2.0D0*FTI(I)/RATIOR)  *PSI  (I) ) 
C 

HPB  (I)  =  -  (CPR  (I)  *  (T  (I)-T  (1*1)  )  +  APR  (I)  *  (T  (I)  -T  (1-1)  )  + 

1C3ARPR  (I)*  (PSI  (I)-PSI  (1  +  1)  )  *  A  B  A  RPR  (I)*  (PSI  (I)  -  PSI  (I- 1)  ) 

2*CQ  (I)  *RO  (I)  *QWW  (I)  *DELZ* (T (I* 1 ) -T  (1-1 ) ) /2.0D0 

3-  (2.  ODO*CTT  ( I)  /RATIO  R)  *T  (I)  ) 

350  CONTINUE 

H  (2)  =  (  (2.0D0*FTT  (2)  /  RATIOS)  -C  (2)  )  *  PSI  (2)  *C  (2)  ♦PSI  (3)  *2. 0  DO  * 
1DELZ*  (EV— KPLUS  (2) )  *C  B  A  R  (2) *  (T  (3) -T  (2) ) 

HPR  (3)  =2.0D0*CTT  (3)  *T(3)/RATIOR  +  (APR  (3)  *CQ  (3)  *RO  (3)  *Q«K  (3)  * 
1DELZ/2-0D0- (ABARPR  (3) *  A  3  A  R  (3)  ) 

2/A  (3)  )  *TI  (2)  »CPH  (3)  *  (T  (4)  -T  (3)  )  +  APR  (3)  *  (T  (2)  -T  (3)  )  +ABARPR  (3)  * 
3ABAR  (3)  /A  (3)  *  (T  (3)  -  T  (2))  +C3ARPR  (3)  *  (PSI  (4)  -PSI  (3))  +2.0D0* 
4DELZ  * A  B A  RPR  (3) /A  (3) *  (EV-K MINUS (3)  ) -CQ (3) *20  (3)  (3)  *DELZ* 

5  (T  (4)  -T  (  2)  )  /2.  ODO 
I=L 

352  H(I)  =H(I)  »C  (I)  *PSI  (IH)  -t-CBAR  (I)  *T(I*  1) 
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HPR  (I)  =HPR  (I)  »CPR  (I)  *T  (1*1)  ♦CBA8PR  (I)  *PSI  (I*1) 

C- — FORA  RHS  :  HPR- (E)  *  (PSI)  =U  ;  H-(F)*(T)=7 
DO  9875  1=2, L 
IF  (I.  EU-  2 )  GO  TO  9874 

V(I)=H(I)-F(I,1)*T(I-1)-P(I,2)*T(I)-F(I,3)*T(I*1) 

GO  TO  9873 

987  4  V  ( 2)  =H  (2)  -  F  (2,2)  *T(2)-F(2,3)  *T  ( 3 ) 

9873  GI  (I) =G  (I,  1) 

GII  (I)=G  (1,2) 

GUI  (I)  =G  (1,3) 

9875  CONTINUE 

CALL  TRIDAG  ( 2  ,  L,  GI ,  GII ,  GUI  ,  V ,  PSII ) 

DO  9870  1=3, L 

U (I) =dPR  (I)  -  E (I,  1)  *  P  SII (1-1) -E  (1,2) *?SII  (I ) -S  (1,3)  *PSII (1+1) 
DI  (I)  =D  ( 1 ,  1) 

DII  (I)  =D  (1,2) 

Dill  (I)  =D  (1,3) 

9870  CONTINUE 

CALL  TRIDAG  ( 3 , L , DI , DII , Dili , U, TI ) 

CALL  THE5P2  (PSII , Ku , THET A I , N6 6 , PSI 66, P66 , L) 
IF(KOUNT.GT.KMAX)  GO  TO  888 
DO  77  1=3, L 

IF  (DABS (TI  (I) -T (I) ) - LE- EPS)  GO  TO  77 
GO  TO  9 
77  CONTINUE 

DO  710  1=2, L 
TO  (I) =T  (I) 

T(I)=TI  (I) 

IF  (I. EQ.  1)  GO  TO  710 
PSIO(I)  =PSI  (I) 

PSI  (I)  =PSII  (I) 

THETAO (I) =THET A (I) 

THETA  (I)  =THETA I  (I) 

CCO(I)  =CC  (I) 

710  CONTINUE 
GO  TO  27  1 

9  DO  711  1  =  2, L 
TO (I) =T  (I) 

T  (I)  =  (TI  (I)  *-TO  (I)  )  /2.0D0 
IF  (I. SQ.  1)  GO  TO  711 
PSIO  (I)  =  PSI  (I) 

PSI  (I)  =  (PSII  (I) +PSIO  (I) ) /2.0D0 
THETAO  (I ) =THET A (I) 

THETA  (I)  =  (THETAI  (I)  » THETAO (I)  )  /2.0DO 
CCO  (I)  =CC  (I) 

711  CONTINUE 
GO  TO  699 

777  DO  17  1=2, L 
CCO (I) =CC (I) 

DPSIO  (I)  =DP SI  (I) 

FO (I )  =  FF  (I) 

LAW  DAO  (I)  =  L  A  N  D  A (I) 

KO  (I )  =  K  (I) 

CLQO  (I) =CLQ  (I) 

QWO  ( I ) =  Q  W  (I) 

DLO  (I)  =DL  (I) 

DTO  (I)  =DT  (I) 

17  CONTINUE 
GO  TO  778 

C - PRINT  OUT  CALCULATED  VALUES 
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271  IP  (K0DE2. EQ.  1)  GO  TO  6016 
IP  (J-KODE4.  NE.  0)  GO  TO  270 
KODE4=KODE2*J 

6016  WBITE (6, 50 19)  J , KOUN T,  (TI (I )  , I  =  2 , L) 

WRITE (6,  5020)  J,  (PSII  (I)  ,1  =  2, L) 

5019  FOHH AT (1  HO, • J= • , 13, 3X, • KOUNT=' ,I3,/(1H0»'T',6G15.7)) 

50  20  FOR  HAT (1  HO, • J  = • , 13,/  (1  HO, ' PS  I* ,6G15.7)  ) 

6015  WRITE(6,700) 

70  0  FORMAT  (1  HO,  •  =  =  =  =  ===  =  ===^==  =  =  =  =  =  =^====  ==  =  =  ========—===  ===  ===' ) 

270  CONTINUE 

WRITE  (6,702)  J 

702  FORMAT {1  HO ,' NORMAL  T ER MI N AT  ION , J= • , 1 5 , 2X , • ITER ATIO NS • ) 

GO  TO  999 

887  WR ITE ( 6 , 87) 

87  FORMAT (1  HO, • WARN ING — PROFILE  DRIED  TO  THE  BOTTOM') 

GO  TO  999 

888  WRITE  (6,889) 

889  FORMAT(1 HQ ,' LIMITS  EXCEEDED;  NO  CON YERG  ENCE* ) 

999  STOP 

END 

— SUBROUTINES 

SUBROUTINE  DTCOEF  (DT V , DT L , D PSI , EV , N R  ,  N I , QCD, QS ,KODE3) 

IMPLICIT  REAL* 8  (A— H,0— Z) 

COMMON/D5S  HE/  PSI  (50) , THETA  (50)  ,T(S0)  ,QW(50)  ,QH(50),DELTAT(50)  , 
1  LAMDA (50 )  , CLQ ( 50 ) , S V P  ( 50 )  , RO  (50 ) , 3 1  ( 50 ) , P 2  (4 , 50) , T 2 ( 50)  , 

2T1 (50)  ,P 1 (4,50) ,T3  (50)  ,P3 (4,50)  ,THET1 1  (50)  , ? 1 1  (4 , 50)  , H V AP ( 50)  , 
3THETA5  (5  0)  ,P5(4,50)  , T9  (50)  , P 9 ( 4 , 50 ) , FF ( 50 )  , 

4PSI 66(50)  ,P66  (4,5  0)  ,T7  (50)  ,P7  (4  ,50)  ,T8  (50)  , 

5 P8  (4,50)  , SR H  (50)  ,K(50)  ,QWO (50)  ,  DAW  (50)  ,DT(S0)  , 

6THETAO  (50)  ,NJ,L,L1,tCL,MW,DELZ 
COMMON/NN/  N1,N2,N3,N5,N7,N8,N9,N11,N66 

DIMENSION  DTV(L)  ,DTL(L),DPSI(L)  , 

1 QC V  (50) ,QCD  (L)  ,QS  (L) 

REAL *8  K ,MW , LAMDA, KS AT, MI,aiTREP,H YAP 
DATA  R/8  3143  2.  0D  +  2/, GR/9 80.  665  DO/ 

DO  333  1  =  2, L 

33  IF  (PSI  (I).GE.O.ODO)  GO  TO  303 

DTL  (I)  =K  (I)  *R*DLOG  (SRH  (I)  )  /  (MW*GR) 

30  DTL  (I)  =D ABS  (DTL(I)  ) 

IF  (KODE3. EQ. 0)  DTL(I)=0.0D0 

DT  (I) =DTL  (I)  +DTV  (I) 

IF  (KODE3 . EQ. 0)  DT(I)=0„0D0 
IP  (I.LE. 3)  GO  TO  331 

QW  (I)  =-DPSI  (I)  *  (  PS  1(1+1)-  PS  I  (1-1))  /  (2.  0D0*DELZ)  -DT  (I)  * 

1  (T  (I  ♦  1 )  -T  (1-1)  )/(2.0D0*DELZ)  (I) 

GO  TO  333 
303  DT (I ) =  0. 0 DO 
DTL  (I)  =DT  (I) 

IP  (I.LE. 3)  GO  TO  331 
IP  (I. EQ. L)  GO  TO  333 

29  8  QW  (I)  =- (2.0  DO/  (1.0D0/K  (1-1)  + 1. ODO/S  (1+ 1) )  ) * ( (PSI (1*1)  -PSI (1-1)  ) 
1/(2.0D0*DELZ)  -  1.0D0) 

GO  TO  333 
331  QW  (2) =EV 
333  CONTINUE 

DO  100  1=2, L 
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QCV  (I)  =  RO  (I)  *CLQ  (I)  *  DELTAT  (I)  *Q»  (I) 

IF  (KODE3. EQ. 0)  QCV(I)=0.0D0 

QH  (I)  =QCD  (I)  ♦QS  (I)  *QCV  (I) 

100  CONTINUE 
440  RETURN 
END 
C 

SUBROUTINE  DPSICF (DT H V  , D PSI , KS AT ,DTH £T A , DP5IV) 

C 

IMPLICIT  RE  A  L*  3  (A— H,0— Z) 

COMMON/D5WHE/  PSI  (50)  , T BETA  (50)  , T  (50) , QS  (50)  , QH  (50 ) , DELTAT (50)  , 
1LAMDA  (50)  , CLQ ( 50 )  , S V P  ( 50 )  , RO (50)  , 3 1 (50 )  , P 2 ( 4 , 50)  , T 2  (50)  , 

2T1  (50)  , P 1 (4 ,50)  ,13  (50)  ,P3 (4, 50)  ,THET1 1  (50)  ,P1 1  (4,50)  ,HVAP(50) , 
3THETA5  (5 0)  , P 5 ( 4 , 50)  , T9  ( 50)  , P9 ( 4 , 50)  , FF  ( 50)  , 

4PSI66(50)  ,P6  6  (4,50)  ,T7  (50)  ,P7(4  ,50)  ,T8  (50)  , 

5P8  (4,5  0)  , SR H  (50 )  ,K  (5  0)  ,  QWO  (50)  , DAW  (50)  , QT(50)  , 

6THETAO  (50)  ,NJ, L,  LI , KL, MW, DELZ 
CO  a HON/S N/  N1,N2,N3,N5,N7,N8,N9,N11,N6o 
C 

DIMENSION  DTHV  (L)  ,DPSI  (L) , KSAT (L) , DT8ET A (L)  ,DPSIV(L) 

REAL *8  K,MW,  LA HD A, KS AT , fl I, MITR EF, H VAP 
DATA  R/3  31  432.  QD ♦ 2/, GR/9 30.  665D0/ 

DO  455  1=2, L 

IF  (PSI  (I)  .GE.0.0D0)  GO  TO  45 
DPSI (I)  =K (I) +  DPS  IV (I) 

IF  (FF(I)  .EQ.0.0D0)  GO  TO  453 
DTHETA  (I)  =DPSI  (I)/FF  (I) 

GO  TO  454 
45  DPSI  (I)  =K  (I) 

453  DTHETA  (I) =0. 0D0 

454  DPSI (I) =DA3S  (DPSI  (I) ) 

455  CONTINUE 
RETURN 
END 

C 

SUBROUTINE  HYCON  D ( J, KSAT,3ITR£F , ROT REF, DTHV, DPSI V  , 
1D£LTRF,DTV,EV) 


IMPLICIT  R  EAL*  8  ( A-H, 0— Z) 

COMMON/D5WHE/  PSI(50) , THETA  (50)  ,T(50)  ,QW(50)  , Q H  (50 ), DELTAT (50)  , 
1LAMDA  (50)  , CLQ (50)  , SVP (50)  , RO (50 )  , MI ( 50 )  , P2 (4 , 50)  , T 2 ( 50)  , 

2T1  (50) ,P  1  (4,50)  , T3 (SO)  ,P3 (4, 50)  , THE! 1 1  (50) ,P1 1  (4 ,5 0)  ,HVAP(50)  , 
3THETA5 (50)  ,P5(4,50)  , T9 (50)  , P9  (4 ,50)  ,  FF  (50)  , 

4PSI66  (50)  ,P6  6  (4,  50)  ,T7  (50)  , P7 (4 ,50) , T3  (50)  , 

5P8  (4,50)  ,  SRH  (5  0)  ,K (50)  , QWO  (50)  , DAW  (50)  ,DT(50)  , 

6THETAO (50) , N J , L , LI , KL , 3 W , DE LZ 
COMMON/NN/  N  1, N2,N3,N5,N7, N8.N9, N1 1,N6  6 

DIMENSION  DTHV  (L)  ,DPSIV  (L)  , 

IDELTRF(L)  ,  K S  AT  (L)  ,  DT  V  (  L)  ,BL(50) 

BEAL* 8  MI,MITREF,K,KSAT,MW, LAMDA,HVAP 
DATA  R/3  31 43  2. OD *  2/, GR/9 80.  665D0 / 

DO  111  1=2, L 

IF  (PSI (I)  .LI.O.ODO)  GO  TO  111 
IF  (DELTRF  (I)  .  G  E.  1.0D0)  GO  TO  102 

101  K  (I) =KSAT (I) 

GO  TO  111 

102  K (I) =KS AT  (I) * (MITREF/ROTREF)  * RO  (I) /MI (I) 

111  CONTINUE 
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IF  (J.LE.  1000)  GO  TO  113 

C - CALCULATION  OP  ONSAT  K  USING  NON-ISO THERMAL  BATES  FLUX  EQN 

DO  112  1=2, LI 

IF  (P3I (I) .LT. 0. 0D0)  GO  TO  103 
GO  TO  112 

103  IF  (I. EQ. 3)  GO  TO  112 

BL  (I)  =B*DLOG  (SBH  (I)  )  / ( MW  *G  B) 

IF  (I.  EQ.  2)  GO  TO  99 

K  (1)  =  (-QH  (I)  -DPSIV  (I)  *  (PSI  (1  +  1)  -PSI  (1-1)  )  /  (2,  0  D0*D  ELZ) 

1-DTV  (I)  *  (T  (1  +  1)  -T  (I-  1)  )/  (2.0D0*DELZ)  ) 

2/(BL (I) * (T (1+1) -T  (I-  1)  )/  (2. 0 DO *DELZ)  ♦ 

3  (PSI (1  +  1)  -  PSI (1-1) )  /  (2.0D0*DELZ) -1 .0D0) 

GO  TO  98 

9  9  K(I)  =  {— QB(I)— DPSIV  (I)*  (PSI  (1  +  1)  — 

IPS  I  (I)  )  /DELZ-DTV  (I)  *  (T  (1  +  1)  -T  (I)  )  /DELZ)/(BL  (I)  *  (T  (1+  1)  -T(I)  )  / 
2DELZ+  (PSI  (1+  1)  -  PSI  (I)  )  /D ELZ- 1.0 DO) 

98  K  (I)  =QAB S  (K  (I)  ) 

1  12  CONTINUE 
113  BETJfiN 
END 
C 

SUBBOUTINE  EVAPO (EV, DELT , DT V , D PSIV , DTL , J , KOUNT , KOD E3  ,  KOD B5) 

C 

IMPLICIT  3EAL*8 (A-H,0— Z) 

CO MM ON/D 5 SHE/  PSI(SO)  , THETA  (50)  ,T(50)  ,QW(50) ,QH(50) , DELTAT (50)  , 
1  LA MD A ( 50 )  ,CLQ  (50)  ,SVP(50)  ,RO  (50)  ,MI  (50)  ,P2  (4,50)  ,T2(50)  , 

2T1  (50) ,P1  (4, 50) ,T3  (5  0)  ,P3  (4, 50)  , THET11  (50)  ,?1 1  (4,50)  ,HVAP(50) , 
3THET A 5 (50)  ,P5(4,S0)  , T9 (50)  , P9  (4 , 50 )  , FF  (50 )  , 

4PSI66  (50)  ,P66  (4,50)  ,T7 (50)  ,P7 (4,50)  ,T8  (50)  , 

5P8  (4,50)  ,SRH  (5  0)  ,K  (50)  ,QWO  (5  0)  , DAH  (50)  ,DT(5Q)  , 

6THETAO (50) , N J , L, L 1 , K L , MB ,DELZ 
COMMON/NN/  Nl,N2,N3,N5,N7,N8,N9,N11,N66 
C 

DIMENSION  DTV  (L)  , DPSIV  (L)  ,DTL(L) 

BE AL*8  MW,K,LAMDA,HVAP,MI 

DATA  3/831432. OD+2/, GH/980. 66500/ 

DT  (3)=DTV  (3)+DA3S(K(3)  *B*DLOG  ( SRH ( 3) ) /  (MM *GR) ) 

IF  (KODE3 . EQ. 0)  DT(3)=0.0D0 

48  QB  (3) =-  (K  (3)  +0PSIV  (3) )  * (PSI (4)  -  PSI (2) ) /(2. 0D0*DELZ) 

1-DT (3) * (T (4)-T  (2) ) / ( 2. 0 D0*DSLZ )  +  K(3) 

IF  (KODE5. GT. 0)  GO  TO  40 

444  EV  =  DABS (QW  (3) ) -  (THETA  (2) -THETAO (2)  +THETA  (3) -THETAO (3)  )  /2-ODO 
1*DELZ/DELT 

IF  (QB  (3)  .  GT.  0 .  0D0)  SV=-Q»  (3)-  (THETA  (2)  -THETAO  (2)  +THETA  (3) 
1-THETAO (3) ) /2. 0D0*D£LZ/D£LT 
45  B RITE ( 6 , 4  7)  J, KOUNT, EV  ,QW (3)  ,K  (3)  , THETA  (2)  , THETA (3) 

4  7  FORMAT (1  HO, • J= ' , 13, 2X,  • KOUNT • , I  3 , 2 X , ' S V= • , G 1 5. 7 , 2X , • QB ( 3)  =  • , 

1G1 5. 7, 2X ,' K  (3) ='  ,G15.7,2X, 'THETA  (2 ,3)  '  ,2G15.7) 

EV=- EV 
GO  TO  41 

40  EV=0. 09D0 

41  RETURN 
END 

C 

SUBROUTINE  HCAPT  ( XM , XO ,  L, THETA , CC , KL) 

IMPLICIT  R  E  A  L*  8  (A— H,0-Z) 

DIMENSION  XM  (KL)  ,XO (KL)  , THETA (KL)  ,CC  (KL) 

DO  213  1=2, L 

CC (I) =0. 46D0*  XM  (I)  +0.60D0*XO(I)  +THETA(I) 

213  CONTINUE 
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RETURN 
EN  0 

SUBROUTIN  E  HFLUX(  DTHV, DTV, QVAP, DPS IV, QS,QCD,KODE3) 

IIIPLICIT  REAL*8 (A-H,0-Z) 

COMMON/D5«HE/  PS  I  (50)  ,  THETA  (50)  ,T  (50)  ,QW  (50)  ,QH  (50)  ,DELTAT  (50)  , 
1  LAND A (50) ,CLQ (50)  ,SVP(50)  , HO  (5  0) ,HI  (50)  ,P2  (4 ,50)  ,T2(50)  , 

2T1  (50)  ,P1  (4,50)  ,1  3  (50)  ,P3(4,50)  , TH  £T 1 1  (50)  ,P11  (4 , 50)  , H V AP  (50) » 
3THETA5 (5 0)  ,  P  5 ( 4 , 50 )  ,  T9 ( 50) , P 9 (4 , 50 ) , FF  (50 )  , 

4PSI6  6  (SO)  ,P6  6  (4,  50)  ,T7  (50)  , P7 (4  ,50) , T8  (50)  , 

5P8  (4,50)  ,SRH  (5  0)  , K  (5  0)  ,QW0(50)  ,  DAW  (50)  ,  DT  (50)  , 

6THETAO (50) ,  NJ, L, LI, KL, MK, DELZ 

COHflON/NN/  Nl,N2,N3,N5,N7,N8,N9,N11,N66 

DIMENSION  DTH  V  (L)  ,  DT  V  (  L)  ,QVAP(L)  ,QCD(L)  ,QS(L)  ,QCV(SO)  , 

1DPSIV (L) 

HE  AL*8  LAUDA, HVAP,KW,K, HI 

DATA  H/831432.  OD + 2/, GR/9 80. 665D0/ 

POB=0. 4 1  DO 
DO  100  1  =  2, L 

IF  (PSI  (I).GE.O.ODO)  GO  TO  99 
Z  ETA= 1 . 8  DO 

IF  (THETA (I) . LE.0.09D0)  GO  TO  95 
RK=  (POR-THETA  (I)  ) / (POR-O . 09D0 ) 

PF=  (2. 0 DO/ 3 . 0 DO)  *  (POR-THETA  (I)  +  RK* THETA (I) ) 

GO  TO  96 

95  PF=2.0D0*POR/3.0D0 

96  RT=KW*GR/ (R*  (T  (I) +273.  16D0) ) 

RTS=RT*RT 

DPSIV  (I)  =  PF*  DA  W  (I)  *  S  VP  (I)  *RTS*S3H  (I) 

97  IF  (KODE3. EQ. 0)  D PS1 V  (I ) = 0. 0 0 0 
GO  TO  98 

99  DPSIV (I) =0. 0D0 
QS  (I) =0. 0D0 

98  IF  (I. EQ. 2)  GO  TO  30 

QCD(I) =  - L AN  D  A  (I) *  (T  (1+  1) -T(I-1)  )/(2.0DO*DELZ) 

QS (I) =-DPSIV  (I)  *  HO (I) +HVAP (I) *  (PSI  (I  +  1)-PSI  (I- 1)  )/ (2.0D0*DELZ) 
SP=  (SVP  (1+  1)  -S  VP  (I-  1)  )  /  (X  (1+  1)  -T  (I-  1)  ) 

GO  TO  40 

30  QCD  (2) =- LAM  DA (I)  *  (T  (1+  1) -T (I) ) /DELZ 

QS(2) =— DPSIV  (I)  *RO (I)  *HVAP(I)  * ( PSI  (I > 1 ) —PSI (I)  ) /DELZ 
SP= (SVP (1+ 1) -S  VP  (I) ) / (T (1+1 ) -T  (I) ) 

40  IF  (PSI  (I  ) -GE.  0.0D0)  GO  TO  331 
IF  (FF  (I)  .EQ. 0. 0D0)  GO  TO  45 
DTHV  (I)  =  DPS IV  (I)  /FF  (I) 

45  IF  (T  (I)  .  GE.  2  5.  OD  0)  ZETA=1.3D0 
56  DTV  (I)  =ZETA*PF*3T*DA  W(I)*SHH(I)*SP 
IF  (KODE3 . EQ. 0)  DTV(I)=O.ODO 
IF (I. EQ. 2)  GO  TO  300 

QVAP  (I)  =-DPSIV  (I)  *  (PSI  (1+  1)  -PSI  (I-  1)  )  /(2.0D0  +  DELZ) 

1-DTV  (I)  *  (T  (1  +  1  )  -T  (I-  1)  )  /  (2. 0  DO*  DELZ) 

GO  TO  100 

300  QVAP  (2)  =- DPSIV  (2)  *  (PSI  (3)  -  PSI  (2)  ) /DELZ  -DTV  (2)  *  (T  ( 3) -T  ( 2)  )  /DELZ 
GO  TO  100 
331  DTV (I) =0 . 0  DO 
DTHV  (I)  =0.  ODO 
QVAP  (I) =0. ODO 
100  CONTINUE 
321  RETURN 
END 
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SUBROUTINE  FCO  EF  (L,PSI , THETA , SS , FP , SI) 

IMPLICIT  R  EA  L«  8 (  A-H  ,  O-  Z) 

DIMENSION  PS  I  (KL) , THETA (KL) , FP (L) 

DO  45  1=2, L 

IF  (PSI (I)  .GE. 0.0D0)  GO  TO  35 
IF  (I. EU- 2)  GO  TO  33 
IF  (PSI  (1*1) .GE.0.0D0)  GO  TO  34 
IF  (PSI  (1*1) . EQ. PSI  (1-1) )  GO  TO  32 

FF  (I)  =  (THETA  (I  *  1 ) -T  H  ST  A  (I- 1 )  )  /  (  PSI  (1*1)  -PSI  (I-  1)  ) 

GO  TO  45 

33  IF  (PSI  (2)  .  EQ.PSI  (3)  )  GO  TO  32 

FF  (I)  =  (THETA  (1*1) -THETA  (I)  )  /  (PSI  (IM)-PSI  (I)  ) 

GO  TO  45 

34  FF  (I)  =  (THETA  (1*1)  -THETA  (I)  )  /  (-PSI  (I)  ) 

GO  TO  45 

32  FF  (I )  =0.  0 DO 
GO  TO  45 

35  PF  (I)  =SS 
45  CONTINUE 

RETURN 

END 

C 

SUBROUTINE  DLCOE F  (DTH 7 , DL ,D PSI V) 

C 

IMPLICIT  REAL*8  (A-H,0— 2) 

C0MM0N/D5HHE/  PSI(50),THETA(50),T(50),QM(50),QH(5Q), DELTA! (50)  , 
1  LAND A  (50)  , CLQ ( 50 ) , S7 P ( 50) , RO  (50)  , MI  ( 50 ) , P2  (4 , 50)  , T2 ( 50)  , 

2T1  (50)  ,P 1 (4 ,50) ,T3  (50)  ,P3  (4,50)  ,THET1 1  (50)  ,P1 1  (4,50)  ,HVAP(50)  , 
3THETA5 (5  0)  , P5(4, 50)  , T9  (50)  , P9 (4 ,50)  , FF (50)  , 

4PSI66  (50)  ,P66  (4,50)  , T7  (50)  , P7  (4  ,50 )  , T3  (50 )  , 

5P8 (4,5  0)  , SRH  (50)  ,K  (50)  ,Q«0  (50) , DAM (50)  ,DT (50)  , 

6THETAO  (50)  ,NJ, L, LI , K L, MB , DELZ 
COMMON/NN/  N1,N2,N3,N5,N7,N8,N9,N11,N66 
C 

DIMENSION  DTHY  (L )  ,DL  (L)  ,DPSIV (L) 

RE AL*8  MM, HV AP, KSAT, K, LAMDA, MI, MITREF 
DATA  R/8  3143  2.  0D*2/, GR/980.  66 5 DO/ 

DO  65  1=2, L 

IF  (PSI  (I).GE.O.ODO)  GO  TO  75 
DL  (I)  =  80  (I)  *H7AP  (I)  «DPSI7(I) 

GO  TO  65 
75  DL (I)  =0- 0D0 
65  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  DENS (T , L , RO , N 1 , T 1 , P 1 , KL) 

IMPLICIT  REAL* 3 (A-H,0-Z) 

DIMENSION  RO  (KL)  , T  ( K  L)  ,  T  1  (N  1 )  ,  PI  (4  ,  N  1 ) 

NR  =  N  1-  1 

CALL  SPLINE  (NR  ,T1, PI ,N 1) 

CALL  CALCCF (NR ,T 1 , P 1 ,N 1) 

DO  63  1=2, L 
X  =  T  (I) 

63  RO  (I) =PCU3IC  (X,NR,T1 ,P1,N1) 

RETURN 

END 

C 

SUBROUTINE  VISCO  (T, L , MI , N2 ,T2 , P2 , KL) 
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IMPLICIT  R£AL*8  (A-H,0-Z) 

DIMENSION  T2  (N  2)  ,P2(4,N2)  ,NI(KL),T(KL) 

RE AL*8  MI 
NR=N2- 1 

CALL  SPLINE  ( N  R  , T  2  ,  P2  ,  N  2) 

CALL  CALCC?  (NR,T2,P2,N2) 

DO  63  1=2, L 
X=T  (I) 

63  MI  (I)=PCU3IC(X,NR,T2,P2,N2) 

RETURN 

END 

SUBROUTINE  V AP SP L (T , L, H V A P , N3 , T3 , P 3 , KL ) 

IMPLICIT  REAL* 8 ( A- H , 0- Z) 

DIMENSION  T  3  (N  3)  ,P3  (4,N3)  ,HVAP  (KL)  ,  T  (K  L) 

RE AL *8  H V AP 
N  R=  N3-  1 

CALL  SPLINE  (NR, T3,P3,N3) 

CALL  CALCC F (NR,T3,P3,M3) 

DO  6  3  I  =  2 , L 
X  =  T  (I) 

63  H  V  AP (I)  =  PCU3IC  (X,NR,T3,P3,N3) 

RETURN 

END 

SUBROUTINE  THE  SP  2  (P SI , KL.THET A , N6 6 , PS  166 , P66 , L) 
IMPLICIT  REAL*8 (A-H.O-Z) 

DIMENSION  THET A ( KL)  , PS  I  (KL)  , P66  (4, N66)  , PSI66 (N  66) 
NR=  N  66- 1 

CALL  SPLINE  (NR,PSI66,P66,N66) 

CALL  CALCCF (NR ,PSI66,P66,N6b) 

DO  64  1=2, L 

IP  (PSI  (I)  .GE.0.0D0)  GO  TO  62 
IF  (PSI  (I)  -  LE.-  1.  0D5)  GO  TO  63 
X=PSI  (I) 

663  THETA  (I)  =  P CUBIC (X,N8,PSI66,P66,N66) 

GO  TO  64 

62  THETA (I) =0. 4  IDO 
GO  TO  64 

63  THETA (I) =0. 037D0 

64  CONTINUE 
RETURN 
END 

SUBROUTINE  S VP SP  L (T , L, S VP , N7 , T7 , P7 , K L) 

IMPLICIT  REAL*8 (A-H,0-Z) 

DIMENSION  T7  (N  7)  ,  P7  (  4  ,  N  7)  ,  S  VP  ( KL)  ,  T  (KL) 

NR=N7- 1 

CALL  SPLINE  (NR ,T7,P7 ,N7) 

CALL  CALCCF(NR,T7,P7,N7) 

DO  63  1=2, KL 
X  =  T  (I) 

63  SVP(I) =?CUBIC ( X,NR,T7,P7,N7) 

RETURN 

END 

SUBROUTINE  CLQCF(T,L,CLQ,N8,T8,P8,KL) 

IMPLICIT  REAL*8 (A-H.O-Z) 

DIMENSION  T  8  (N  8)  ,Pd(4,NS)  ,CLQ  (KL)  ,T  (KL) 

NR-N8- 1 
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CALL  SPLINE  (NR  ,T8,P8 ,N8) 

CALL  CALCCF(NR,T8,P8,N8) 

DO  63  1=2, L 
X=  T  ( I) 

63  OLQ (I) =  P  CUBIC ( X, NR,T8,P8,N8) 

RET JRN 
END 

SUBROUTINE  S PLI N E  ( N R , X  I , P , N S IZ E) 

IH1PLICIT  R  E  A  L  *  3  (  A  -  H ,  0-  Z ) 

DIMENSION  XI  (NSIZE)  ,  P  (  4,  NSIZE)  ,  D  (5  0 )  ,  D I  AG  ( 50) 

DI AG (1 ) = 1. 0D0 
D  ( 1)  =0 . 0  DO 
NP1=NR+ 1 
DO  10  M=2,NP1 
D  (H) =XI  (M) -XI (M-  1) 

10  DIAG  (M)  =  (P  (1  ,M)-P  (1,«-1)  )/D(H) 

DO  20  M=  2 ,  N  R 

P (2,M) =3. 0  DO*  (D  CM) *DIAG  (M  +  1) +  D (M+1 ) *DIAG (B)  ) 

20  DIAG  (M)  =2.0D0«  (D  (M)  +  D  (M+1)  ) 

DO  30  M=  2 , NR 
G  =  -D  (fl  +  1 )  /DIAG  («-  1) 

DIAG  (M)  =  DIAG  (M  )  +-G*D  (M-  1) 

30  P  (2,M)=P  (2,M)+G*P(2,M-1) 

NJ=NP1 

DO  40  M=  2, NR 
NJ  =  NJ-  1 

4  0  P  (2,NJ) =  (P  (2,NJ) -D(NJ) *P  (2, NJ+1)  ) /DIAG (NJ) 
RETURN- 
END 

SUBROUTINE  CALCCF  (NR, XI, P, NSIZE) 

IMPLICIT  REAL*3  (A-H,0-Z) 

DIMENSION  XI  (NSIZE) , P (4, NSIZE) 

DO  10  1=  1  , NR 
DX=XI  (I*  1 )  -XI  ( I) 

DIVD  F  1  =  (P(1,I+1)-P(1,I))  /DX 
DIVDF3=P  (2,1) +  P (2, I* 1) —  2.0D0*DIVDF 1 
P  (3,1)  =  (DIVDF1-P  (  2 , I) —  DI VDF3) /DX 
P  (4,1) =DI VDF3/DX/DX 
10  CONTINUE 
RETURN 
END 

FUNCTION  PCUBIC (XBAR , NR , XI , P, NSIZE) 

IMPLICIT  REAL*8 (A-H,0-Z) 

DIMENSION  XI  (NSIZE)  , P  (4 , NSIZE) 

1=  1 

DX=XBAR-XI  (I) 

IF  (DX)  10,30,20 
10  IF  (I.EQ.  1)  GO  TO  30 
1=1-1 

DX  =XBAR- XI  (I) 

IF  (DX)  10,30,30 

19  1=1+1 
CX  =  D  D  X 

20  IF  (I.EQ. NR)  GO  TO  30 
DDX=X3AR-XI (1+ 1) 

IF  (DDX)  30 ,  19  ,  1  9 

30  P CU BIC=  P  (1  ,1)  +  DX  *  (P  (2,1)  +  DX  *  (P  (3,1)  +  DX  *  P (4,1)  )  ) 
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RETURN 
EN  D 

SUBROUTINE  THIDAG(I1,L,AZ,BZ#CZ,DR,VZ) 

IMPLICIT  REAL»3  (  A-H, O-Z) 

DIMENSION  AZ  (50)  ,BZ(50),CZ(S0)  ,D8(50),VZ(50)  ,3  ETA (S0),GAMMA(50) 
BETA  (II)  =BZ  (I  1 ) 

GAMMA  (II)  =03  (I  1)  /BETA  (II) 

I1P1=I1+1 
DO  1  1  =  1  IP  1  ,  L 

BETA  (I)  =  8Z  1 1 )  -  AZ  (I)  *CZ  (1-1)  /BETA  {1-1) 

1  GAMMA  (I)  =  (DR  (I) -AZ (I) *  GAMMA  (1-1 )  ) /BETA  (I) 

VZ  (L)  -GAMMA  (L) 

LAST  =  L-I  1 
DO  2  K  = 1 , LAST 
I=L-K 

VZ  (I)  =  GAMMA  (I)-CZ(I)  * VZ  (I* 1 ) /BETA (I) 

2  CONTINUE 
RETURN 
END 

SUBROUTINE  THCON  D(L,  LA MDA, PSI, THETA , KL , TH ET 1 1 , P 1 1,  Nil) 

IMPLICIT  S  EA  L*  8 ( A-H , O- Z) 

DIMENSION  LAM DA (KL) , PSI (KL)  , THETA  (KL)  ,  THET1 1 (N 1 1)  , PI  1 (4, N1 1) 
REAL*8  LAM  DA 
NR=N 1 1-1 

CALL  SPLINE  (NR ,THET  1  1 , PI 1,N 1 1) 

CALL  CALCCF (NR,THET11 ,  P 1 1 ,  N  1 1) 

DO  63  1=2, L 
X=  TH  ET A (I) 

63  LAMDA (X) =PCUBIC(X,NR,THET11,P1 1 , N 1  1) 

RETURN 

END 

SUBROUTINE  D AWC?  (T , L , D A U , N 9 , T9,  P9 , KL) 

IMPLICIT  REAL*3  (  A-H,  O-Z) 

DIMENSION  T9  (N  9)  ,P9(4,N9)  ,  DA  W  (KL)  ,  T  ( KL) 

N  R  =  N  9- 1 

CALL  SPLINE (NR , T9, P9 ,N9) 

CALL  CALCC? (NR,T9,P9 ,N9) 

DO  63  1=2, L 
X=T  (I) 

63  DAM (I) =PCUBIC(X,NR,T9,P9,N9) 

RETURN 

END 

SUBROUTINE  KSP LE N  (TH ET A , K, N5 ,THET A5 , P5 , L, PSI, T , MI, RO , 

1T1 ,T2,P1 , P2 , N 1 ,N2,KL,KSAT) 

IMPLICIT  REAL*8 (A-H, O-Z) 

DIMENSION  THETAS  (N  5)  ,  P  5  (4  ,  N5)  ,  THETA  (KL)  ,K  (KL)  ,  PSI  (  KL)  ,  MI  ( K  i_)  , 

1  T  1  (N  1)  ,  T  2  (  N  2 )  ,  PI  (  4  ,  N  1 )  ,  P2  (4,  N2)  ,T(KL)  ,  RO(KL)  ,  KS  AT  (KL) 

REAL*8  K , M I, K3 AT 
N  8=  N  5-  1 

CALL  DENS  (T,L,30,N1,T1 ,P1 ,KL) 

CALL  VISCO  (T , L,M I ,N2 ,T2, P2, KL) 

59  CALL  SPLINE  (NR  , THETA  5, P5,N5) 

CALL  CALCCF(NR, THETAS, P5,N5) 

DO  63  1=2, L 

62  IF  (PSI  (I)  .GE.O.ODO)  GO  TO  64 
X  =  TH ET A  (I) 


K  (I) = PCD DIC  (X, NR, THETA  5, ?5,  N 5) 

C -  TEMPERATURE  ADJUSTMENT  OF  K-UNSAT  BASED  ON  TREF=20C 

60  IF  (DABS  (T  (I) -  20.  ODD)  .  G  E.  1.0D0)  K  (I )  =  K  ( I )  *  36.  1 3  6700  *R0  ( I)  /M I  (I) 
GO  TO  63 

64  K  (1) =K5AT  (I) 

63  CONTINUE 

65  RETURN 
END 
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E5  Sample  Injout  Data  and  Output  Results 


A  sample  input  data  to  the  present 
follows;  also  a  sample  output  data  using 


simulation  program 
these  input  data- 


Input  data  foe  the  simulation  program.  The  line  numbers 
below  refer  to  the  Program  Listing  of  Section  E4* 
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2  KSAT(I) =CSAT2 
NN  2=  NI ♦  1 

DO  3  1=  N  N  2 , L  1 

3  KSAT(I) =CSAT3 

4  KSAT  (L) =CS AT4 

CALL  THE3P2  (P SI , K L, TH ETA , N 66 , PSI66 ,P6 6 , L) 

C - 

DO  6  1=2, L 

IF  (I.GT. BO)  GO  TO  7 

RZAD(5,1111)  X B  (I)  , XO ( I) 

GO  TO  6 
7  XM (I) =xa (BO) 

XO  (I)  =XO  (HO) 

6  CONTINUE 

1111  FOSHAT  ( F  1 0. 3 , F 1 0. 5) 

BEAD (5,2  222)  (T9  (I ) , P9  ( 1 , 1)  ,  1= 1 , N9 )  , P9  (2 , 1 ) , P9  ( 2 , N 9) 

BEAD (5, 2 2 22)  (T8  (I)  ,  P8  ( 1 , 1)  , 1= 1 , NS )  ,  P3  (2,  1)  ,P8  (2, NS) 

READ  (5,2  2  22)  (T7  (I)  ,P7  (1,1)  ,1  =  1,  N7)  ,P7  (2,  1)  ,P7  (2,N7) 

BEAD  (5,2222)  (T3  (I)  ,  P3  ( 1 , 1)  ,  1=  1  ,  N3)  ,  P3  (2,  1 )  ,  P 3  (2,  M  3) 

2222  FORMAT  (2F10.  5) 

READ (5, 2  222)  (T2  (I)  ,P2  (1,1)  ,I=1,N2)  ,P2  (2, 1),P2  (2,N2) 

RE AD (5, 2 2 25)  (THETA 5  (I)  , P5  ( 1 , I )  , 1= 1 , N5) , P 5  (2 , 1 )  , P5  (2 , N5) 
2225  FORMAT (F10. 5, DIO. 6) 

DO  250  1=2, L 

Z=I-2 

Z=DELZ*Z 

211  i RITE  (6,260)  I , Z , PSI  (I)  , THETA (I)  ,T (I) 

250  CONTINUE 

260  FOB 8 AT  (120, 6X,F20.8, 4X  ,3F20. 8) 

T(KL)=T(L) 

PSI  (KL)  =  PSI  (L)  -DE-L2*  (1  . 0  DO*  i/KS  AT  (L)  ) 

KODE4=KO  DE  2 

C - BEGIN  TIME  ITERATIONS 

40  DO  270  J= 1 , N J 

IF  (NJ.GT. NSTOP)  GO  TO  999 
KOUNT=0 

IF  (J .  EQ.  1)  GO  TO  693  1 
GO  TO  699 
6981  DO  697  1=2, L 

THETAO (I ) =THET  A (I) 

TO  (I)  =T  (I) 

PSIO (I) =PSI  (I) 

697  CONTINUE 

C - CALCULATE  TRANSPORT  COEFFICIENTS 

699  CALL  HCAPT (XM,XO,L, THETA, CC.KL) 

IF  (J. GE. 43)  GO  TO  999 

CALL  THESP2  (P SI , KL, TH ET A , N6 6 , PSI66 ,P6 6 , L) 

DO  5325  1=2, L 
IF(J.EQ.I)  CCO(I)=CC(I) 

DELTAT  (I) =DABS  (T  (1+1  )-T  (I) ) 

CTT  (I)  =  (CC (I) >CCO (I)  ) /2.0D0 
DELT3F  (I)=DA8S  (T  (I)  -TREE) 

IF  (KODE3. EQ. 0)  DL(I)=0.0D0 

IF  (KODE3. SQ. 0)  CLQ(I)=0.0D0 

359  IF  IPSI (I) -GE. 0.0D0)  GO  TO  361 
EN  =  MW*GR*PSI  (I) 

ED=  R* (T (I) ♦  2  73  .  16D0) 

8E=EN/ED 

SBH (I) =DEXP  (RE) 

GO  TO  5325 
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